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EVALUATION 


For  the  past  few  years,  RADC  lias  concentrated  its  activities  on  critical 
problem  areas  within  the  field  of  analytical  photogrammetry . In  order  to 
maintain  the  technology  base  which  has  successfully  provided  solutions  to  the 
above,  a continual  reexamination  of  the  state-of-the-art  is  essential.  Under 
a modest  level  of  effort,  DBA  Systems  Inc.  lias  developed  a snapshot  of  current 
capabilities  and  has  provided  a roadmap  of  developmental  activity.  Related 
technology  of  interest  to  the  mapping  community  such  as  the  Global  positioning 
system,  improved  inertial  navigation  systems  and  others  have  also  been  drawn 
upon  for  improving  current  triangulation  procedures.  Recommendations  offered 
by  the  contractor  are  being  seriously  entertained  for  pursuit  under  RADC's 
exploratory  research  program  and  in  consonance  with  TPO  2. 


/'  JOHN  R.  CALLANDER 
Project  Engineer 
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1.0  INTRODUCTION 

The  objective  of  this  study  is  to  identify  new 
techniques  in  analytical  photogrammetry  which  offer  significant 
improvements  in  target  location  accuracy  and/or  drastically  reduce 
data  reduction  time  for  lame  data  cases.  In  order  to  satisfy 
these  broad  objectives, DBA  has  reviewed  the  current  analytical 
triangulation  procedures  in  use  by  the  USAF,  DMA  and  other 
Governmental  and  Commercial  organizations.  mhe  results  of  this 
review  are  presented  in  Section  2 in  the  form  of  a basic  analytical 
triangulation  model  directed  at  the  use  of  Aerial  Photography. 

Once  the  current  procedures  were  broadly  defined, 

DBA  proceeded  to  analyze  each  step  to  determine  the  major  short- 
falls and  attempts  to  constructively  identify  alternative  tech- 
niques to  alleviate  these  shortfalls.  The  results  of  this  analysis 
are  presented  in  Section  3. 

Sections  2 and  3 provide  background  and  the  basic 
organizational  framework  for  Section  4 which  discusses  in  detail 
the  alternative  techniques.  Some  of  the  techniques  described  in 
Section  4 are  very  detailed,  whereas  others  are  only  conceptual 
and  would  require  more  in-depth  studv  before  their  ultimate 
utility  can  be  evaluated. 

The  conclusions  and  recommendations  resulting  from 
this  study  are  presented  in  Section  5. 

The  primary  shortfall  which  is  identified  in 
Section  3 relates  to  the  inability  of  most  aerotriangulation 
programs  to  adequately  account  for  systematic  errors  in  the 
observations.  Specific  examples  of  these  errors  are  film  defor- 
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mation,  camera  calibration  parameters,  and  most  importantly, 
systematic  errors  in  auxiliary  observations  such  as  GPS  position- 
ing, inertial  survey  data  and  precise  attitude  sensor  data  etc. 

There  are  a number  of  alternative  techniques  available  for  modeling 
these  systematic  errors.  These  alternative  techniques  are  discussed 
in  Section  4 of  this  report.  Some  of  these  techniques  create 
computational  problems  in  that  new  algorithms  are  required  to 
efficiently  solve  the  least  squares  problem.  Of  particular 
interest  is  an  efficient  algorithm  to  solve  banded-bordered  normal 
equations  which  have  large  highly  patterned  borders. 

Other  shortfalls  which  are  considered,  relate  primarily 
to  increasing  the  overall  throughput  of  the  entire  point  positioning 
process,  not  simply  reducing  the  data  reduction  time.  Included  in 
this  area  are  techniques  for  identifying  available  coverage  of  an 
area,  automatic  and  semi-automatic  editing  schemes,  and  automatic 
data  reordering  algorithms. 

Another  broad  area  which  is  considered  in  this  report 
is  the  incorporation  of  object  space  constraints  in  a aeneral 
photogrammetric  reduction.  This  discussion  is  not  limited  to  the 
normal  distance,  azimuth  and  equal  elevation  constraints  (although 
they  are  discussed) , but  includes  a general  technique  for  including 
geometric  type  constraints  with  specific  reference  to  digital  data 
collection  schemes. 

Naturally  all  items  discussed  do  not  receive  equally 
detailed  development,  however,  sufficient  information  is  presented 
such  that  the  general  approach  can  be  understood  for  all  areas. 

In  some  cases  detailed  mathematical  developments  are  presented, 
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whereas,  in  others  only  a conceptual  analysis  is  presented. 

Those  areas  discussed  in  conceptual  terms  fall  into  two  broad 
categories.  They  are  either  standard  techniques  which  have 
been  documented  elsewhere  and  not  implemented  in  an  aerotrianqu- 
lation  reduction,  or,  the  concept  is  very  new  and  there  wer^ 
insufficient  resources  to  develop  the  details  completely  for 
presentation  at  this  time. 
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2.0 


CURRENT  ANALYTICAL  AEROTR I ANGULATION 


For  the  purpose  of  this  report,  we  will  define  those 
techniques  for  precisely  determining  the  position  of  targets 
identified  on  aerial  photographs  as  an  Analytical  Aerotrianqulation 
System.  Figure  2.0-1  diagrammatically  illustrates  this  system 
and  identifies  ten  major  processes  of  the  system.  These  ten 
nrocesses  can  further  be  broken  into  a number  of  steps,  each  of 
which  can  be  accomplished  in  a variety  of  ways.  The  techniques 
employed  during  any  one  step  often  have  a direct  influence  on 
the  techniques  to  be  employed  during  other  steps  or  subseouent 
processes . 

The  interelationship  between  the  processes,  steps 
and  techniques  of  this  system  operates  in  both  directions.  It 
is  obvious  that  alternatives  selected  at  each  step/process  can 
define  or  limit  the  alternatives  available  for  subsequent  steps/ 
processes.  Experience  has  dictated  that,  although  not  as  obvious 
and  frequently  overlooked,  the  alternatives  selected  for  sub- 
sequent processes  may  also  influence  selection  of  alternatives 
for  prior  processes.  That  is  to  say,  a prior  knowledge  of  the 
constraints  on  a final  product  may  define  the  methods  by  which 
the  product  is  attained.  This  is  particularly  true  for  economical 
and  timely  satisfaction  of  requirements  in  target  location. 

As  a basic  example  of  this  reverse  relationship, 
let  us  assume  that  accuracy  specifications  for  point  locations 
have  been  defined.  In  order  to  meet  these  specifications  one 
must  consider  the  block  adjustment  technique  to  be  employed  as 
well  as  the  scale  of  the  photography  required  to  accommodate 
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FIGURE  2.0-1 

■QGRAMMETPIC  target  location  process 


that  particular  technique.  In  referrinq  to  Figure  2.0-1,  we 
see  that  the  block  adjustment  process  falls  late  in  the  system 
flow,  however,  the  adjustment  parameters,  along  with  scale, 
must  be  known  during  the  coverage  planning  process  such  that 
adequate  photography  will  support  the  adjustment. 

The  following  subsections  present  brief  descriptions 
of  each  process  of  comprising  the  Analytical  Aerotriangulation 
System. 

2 . 1 Area  Selection 

In  most  instances  the  selection  of  an  area  for  target 
location  is  based  on  military  or  political  decisions  and  not  on 
photogrammetric  factors.  It  is  the  science  of  photogrammetry 
which  supports  target  area  selection  and  not  target  area  selection 
that  supports  photogrammetry. 

Photogrammetric  factors  do,  however,  play  a very 
important  role  in  area  selection  if  one  is  planning  to  calibrate 
or  evaluate  a photogrammetric  system  or  technique.  In  this 
circumstance,  certain  parameters  must  be  defined  and  conditions 
satisfied.  For  calibration,  such  parameters  may  consist  of 
adequate  ground  control  and  the  aeometry  of  the  ground  control . 

For  evaluation,  in  addition  to  ground  control,  one  must  consider 
various  altitudes,  terrain  elevations,  slopes  and  auxiliary  data 
acquisition.  Such  considerations,  in  addition  to  others,  are 
discussed  in  detail  in  a previous  report  (DBA,  1974). 

The  Area  Selection  process  is  not  analyzed  in  this 
report  and  is  presented  merely  as  a logical  introduction  to  the 
total  Aerotriangulation  System. 


2.2 


Material  Availability 


Upon  selection  of  an  area  of  interest,  one  must 
determine  what  support  material  is  available  and  must  analyze 
the  suitability  of  any  such  material.  Currently, material 
availability  is  determined  through  incomplete  files  of  oast 
flights  and  the  experience  and  knowledge  of  personnel  assigned 
to  make  this  determination.  Analysis  is  based  on  scale,  coverage 
clarity  (resolution),  amount  of  cloud  cover,  completeness  of 
coverage  ("Holiday"  areas),  etc.  The  majority  of  this  process 
is  dependent  on  human  judgement  and  reliability. 

2 . 3 Coverage  Planning 

Should  the  material  availability  process  prove 
that  coverage  is  inadequate,  the  process  of  coverage  planning 
is  initiated.  The  following  factors  must  be  considered  during 
this  step: 

1.  Scale  of  photography 

2.  Ground  resolution 

3.  Photo  Sensor  type  and  availability 

4.  Control  availability 

5.  Auxiliary  data  collection  system  tvne  and 

availability 

6.  Climatic  considerations 

7.  Time  considerations 

Accuracy  requirements  for  taraet  location  freciuently  dictate 
the  scale  of  the  photography , and  to  some  extent  the  resolution 
specifications.  The  planning  of  missions  must  include  an  analysi 
of  what  types  of  sensors  are  available  for  acquiring  the  data. 
Included  in  the  sensor  availability  must  be  considerations  for 


the  logistics  of  using  then.  The  sensor  availability  process  is 
often  limited  to  a single  sensor,  as  it  may  be  the  only  sensor 
which  can  collect  the  required  data  over  the  specific  area.  This 
default  decision  then  leads  to  the  development  of  techniques 
for  exploitation  of  a specific  sensor. 

Also  included  in  this  evaluation  must  be  a consider- 
ation for  the  availability  of  a data  reduction  program  which  can 
accommodate  this  sensor  and  obtain  the  required  accuracy.  Add- 
itional thought  must  be  given  to  the  calibration  of  the  sensor 
and/or  its  calibratability  during  the  coverage  planning  process. 

Another  area  which  must  be  considered  durincj  this 
planning  process  is  the  availability  of  auxiliary  exposure  station 
data.  As  with  sensor  availability,  consideration  must  be  given 
to  the  ability  to  exploit  this  data  once  it  has  been  gathered. 

The  available  ground  control  must  also  be  considered 
in  planning  the  acquisition  of  photographic  coverage.  This 
consideration  must  include  the  type,  accuracy  and  distribution 
of  such  control. 

Climatic  considerations  are  of  primary  concern  in 
planning  photographic  missions.  Although,  there  is  no  control 
over  this,  the  mission  can  be  scheduled  for  the  most  promisinq 
time  periods  and  to  accommodate  reasonable  contingencies. 

Although  listed  last,  probably  the  most  important 
factor  in  planning  coverage  missions  is  the  time  factor  (i.e. 
when  is  the  requirement  to  be  satisfied).  This  can  have  a 
direct  bearing  on  the  sensor  used,  on  the  control  used  and 
other  factors. 
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As  can  be 
planning  process,  the 
understanding  of  them 
and  an  indirect 


seen  from  this  brief 
techniques  available 
have  a direct  effect 
subsequent 


influence  on  all 


review  of  the  coverage 
and  the  planners' 
on  the  decisions  made 
processes . 


2 . 4 Acquisition 

Material  acquisition  pertains  to  the  required  photography 
as  well  as  to  the  auxiliary  data  and  around  control.  For  the  purpose 
of  this  report,  film  processing,  will  be  included.  An  initial 
review  of  the  coverage  is  nerfomed  to  determine  utility  of  the 
film  for  the  specified  project.  During  this  process,  all  decisions 
made  during  coverage  planning  are  combined  into  a single  system. 

This  timely  combination  of  eauipment , personnel  and  techniques 
may  be  dictated  by  weather  conditions. 


2.5  Point  Selection,  Identification  and  Mensuration 

Once  the  photoqraphic  imagery  is  available,  the  next 
step  in  the  photogrammetric  targeting  problem  is  the  selection, 
identification  and  mensuration  of  the  control  points,  nass  points 
and  target  points. 

The  ground  control  must  be  positively  identified  on 
the  photog-aph.  Typically  this  is  a manual  operation  which 
requires  the  operator  to  stereoscopically  view  a pair  of  over- 
lapping photographs  and  either  measure  the  point  stereoscopical ly 
or  mark  the  imagery  for  subsequent  measurement. 

An  alternative  used  occasionally  is  the  siqnali zation 
of  the  ground  control.  Some  experiments  have  also  been  performed 
with  signalized  pass  points  (i.e.  non-control  points).  In 
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general,  this  is  impractical  due  to  the  time  renuirements 
and  inaccessibility  of  the  area. 

2 . 6 Data  Preprocessing 

Before  a least  squares  adjustment  can  be  initiated, 
the  data  must  be  suitably  transformed,  generated,  edited  and 
arranged  so  that  efficient  algorithms  can  be  utilized  for  the 
least  squares  process.  The  specific  techniques  and  models 
utilized  in  this  phase  must,  of  necessity,  be  matched  to  the 
adjustment  algorithm.  For  this  reason,  computational  routines 
for  data  preprocessing  are  frequently  appended  to  the  "front 
end"  of  the  least  squares  adjustment  program.  However,  for 
economic  reasons  , data  preprocessing  algorithms  are  just  as 
often  small  programs  which  are  executed  independent  of  the 
adjustment.  (It  is  often  more  easy  to  process  a number  of  small 
computer  runs  through  a central  computer  facility  than  a long 
run  and  it  is  easier  to  edit  data  and  rerun  a small  program  than 
a large  one) . 

In  order  to  fully  understand  the  data  preprocessing 
phase  of  an  Aerotriangulation  one  must  first  understand  the 
data,  the  mathematical  model  and  the  least  squares 
algorithm  implemented  in  the  form  of  a block  adjustment  proqram. 

There  are  four  basic  steps  required  in  the  preprocessing 
phase  of  an  adjustment.  These  are: 


1 . 

Data 

transformation , 

2. 

Generation  of 

initial  approximations 

3. 

Data 

editing, 

and 

4. 

Data 

arranqome 

nt. 
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The  data  transformation  steo  converts  the  observations 


and  associated  weiqht  matrices  to  a rc  ference  svstem  con- 
sistent with  the  mathematical  model  utilized  bv  the  block  adjust- 
ment . 

The  generation  of  initial  approximations  for  all 
parameters  (unknowns)  is  required  because  in  general  the  photo- 
grammetric  observation  equations  are  non-linear  with  respect  to 
the  parameters  and  must  be  iteratively  solved  using  a linearized 
set  of  equations.  The  basic  equations  are  linearized  about  some 
arbitrary  initial  approximations.  Corrections  to  the  initial 
approximations  are  computed  and  applied  to  these  arbitrary  values 
and  the  solution  respectively  performed  until  the  corrections 
become  sufficiently  small.  Obviously,  the  closer  the  values  of 
these  arbitrary  initial  approximations  are  to  the  true  value  the 
fewer  the  number  of  iterations  will  be  required.  The  problem  as 
far  as  data  preprocessing  is  concerned,  then,  is  one  of  obtaining 
the  "best"  set  of  initial  values  as  possible  for  all  parameters. 

The  obvious  solution  is  to  set  the  initial  value  equal  to  the 
observed  value  wherever  possible,  however,  for  some  classes  of 
parameters  it  is  more  difficult  to  come  up  with  observations  than 
it  is  to  simply  "guess"  at  initial  values.  Initial  approximations 
for  still  other  classes  of  parameters  can  be  easily  computed  based 
upon  the  use  of  approximate  mathematical  models. 

The  detection  and  elimination  of  aberrant  observations 
(data  editing)  is  performed  durinq  this  portion  of  aerotr iangulat ion 
to  minimize  the  run  time  (and  rerun  time)  for  the  actual  block 
adjustment  program.  This  data  editing  is  frequently  hampered  by 
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the  quality  of  the  initial  approximation  and  the  fact  that  some 
initial  approximations  are  actually  derived  from  the  observations. 

The  last  step  in  data  preprocessinq  is  the  arrangement 
and  storage  of  all  the  data  in  a format  and  on  a media  suitable 
for  efficient  solution  of  the  adjustment.  To  fully  comprehend 
this  step  one  must  be  intimately  familiar  not  onlv  with  the 
mathematical  model  and  the  least  squares  process,  but  also  the 
specific  implementation  algorithms  included  in  the  computational 
reduction.  For  discussion  herein  it  is  assumed  that  one  is 
familiar  with  the  least  squares  solution  and  the  specific  nature 
of  normal  equations  and  the  manner  in  which  they  can  be  reduced  to 
banded  structure  (Gyer  et.al,  1969) . The  problem  of  arranging 
the  data  becomes  one  of  logically  orderinq  the  photographs  in 
a block  to  minimize  the  size  of  the  bandwith  associated  with 


this  reduced  set  of  normal  equations. 


2.6.1  Reordering 

An  analytical  aerotr iangulation  is  based  on  the 
solution  of  a system  of  linearized  normal  equations,  the  formation 
of  which  will  bo  discussed  later  in  this  report.  The  parameters 
to  be  solved  are  usually  iterative  corrections  or  adjustments  to 
orientation  parameters  of  the  individual  photoqraphs.  In  general, 
the  coefficient  matrix  for  this  svstem  is  diagonally  banded  and 
symmetrical.  The  elements  along  the  main  diagonal  are  coefficients 
of  the  corrections  to  the  orientation  parameters  and  the  off  diagonal 
elements  represent  correlations  between  the  parameters. 

It  lias  been  found  that  the  speed  and  efficiency  with 
which  these  normal  equations  can  be  solved  is  directly  related  to 
the  bandwidth  of  the  coefficient  matrix;  the  bandwidth  being  the 
furt herest  distance  of  a correlation  element  from  the  main  diagonal. 

When  the  orientation  parameters  are  grouped  by  individual  photo, 
as  they  usually  are,  the  correlation  elements  furthest  from  the 
diagonal  represent  correlations  between  parameters  of  different 
photos.  This  correlation  is  a consequence  of  the  photos  having 
imaged  a common  object  point.  In  order  to  minimize  the  bandwidth 
and  maximize  the  speed  of  the  solution,  it  is  desirable  to  arrange 
the  parameters  in  the  matrix  so  that  those  for  highly  correlated 
photos  (having  common  image  points)  are  as  close  together  as 

possible.  The  initial  arrangement  of  the  parameters  is  based  on  • j 

a pre-established  ordering  or  numbering  of  the  photos.  If  this 
does  not  provide  the  smallest  possible  bandwidth,  some  rearranging 
or  reordering  is  desirable. 

The  initial  ordering  of  the  photos  is  normally 
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MH 


established  by  a sequential  ohoto  numbering  scheme.'  introduced 


during  mensuration.  The  scheme  used  is  usually  devised  for  con- 
venience in  mensuration  and  is  not  always  the  optimum  for  the 
solution.  The  most  common  scheme  is  down  strin  ohoto  numbering 
in  which  the  photos  in  each  strip  are  numbered  consecutively  in 
the  same  direction  and  the  numbering  sequence  is  continued  from 
one  strip  to  the  next. 

Down  strip  numbering  has  oroven  to  be  ontimum  only 
for  verv  rociular  rectangular  blocks  in  which  the  number  of  photos 
per  strip  (p)  and  the  number  of  strip(s)  are  related  by  ps2s-l. 

For  shorter  broader  blocks  in  which  p-2s-l,  reordering  into  a 
cross  strip  numbering  scheme  is  desirable.  When  more  complex 
correlations  are  introduced  through  cross-strips  or  other 
irregularities,  reordering  by  more  intricate  techniques  is 
required . 

Most  commercial  and  Government  block  adjustment  programs 
require  an  intuitive  assessment  of  the  best  numberina  scheme  and 
manual  ordering  or  reordering  of  the  input  data.  Some  automation 
and  sophistication  has  been  introduced  in  a few  instances.  SAPGO 
(Wong  & Elphingstone , 1972)  has  incorporated  an  automatic  cross- 
strip numbering  capability.  In  order  to  accommodate  some  irreg- 
ularities between  strips,  cross-strip  numbering  is  accomplished 
by  setting  up  a series  of  equally  spaced  parallel  lines  per- 
pendicular to  the  flight  line  of  the  first  strip.  The  parallel 
lines  are  spaced  according  to  the  average  base  length  of  the  first 
strip  and  are  used  as  numberina  zones.  All  of  the  exposure  stations 
that  fall  within  a zone  are  numbered  in  a seguential  manner.  An 
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even  more  universal  technique  has  been  incorporated  into  proaram 
COMBAT  (Brown,  1974).  This  technique  analyzes  the  off-diagonal 
correlation  elements  in  the  actual  coefficient  matrix.  In 
minimizinq  the  bandwidth,  the  longest  rows  in  the  matrix  are 
moved  to  the  center  of  the  matrix  and  the  shortest  rows  to  the 
edge.  The  algorithm  checks  subsequent  rows  to  see  if  they  should 
be  interchanged.  An  interchange  will  take  place  only  if  it  will 
move  a longer  row  closer  to  the  center  of  the  matrix.  With  a 
matrix  of  variable  bandwidth  the  final  reordered  matrix  would 
emerge  appearing  in  a diamond  form.  Experience  shows  that  this 
technique  consistently  produces  the  minimum  bandwidth  ordering 


scheme . 


2 . 7 Block  Adjustments 

There  are  four  basic  categories  of  Analytical 
Aerotr ianqulation  as  represented  in  Table  2.7-1.  The  first 
class  is  intended  to  encompass  all  methods  that  entail  piece- 
wise  (as  opposed  to  simultaneous)  processes  of  adjustment.  It 
includes  most  methods  for  computational  connection  of  independ- 
ently established  models  (i.e.,  semi-analytical  methods)  as  well 
as  methods  involving  numerical  relative  orientation  (i.e., 

’ cant i levering ' ) of  successive  photos  within  each  strip  followed 
by  reconciliation  of  discrepancies  within  and  between  strips  by 
various  polynomial  transformations. 

Methods  of  Class  2 involve  the  assembly  of  independently 
established  models  by  means  of  the  simultaneous,  least  squares 
determination  of  the  particular  sets  of  similarity  transformations 
that  best  interconnect  overlapping  models  while  at  the  same  time 
providing  the  best  overall  accommodation  of  the  block  as  a whole 
to  available  control.  (A  similarity  transformation  performs  trans- 
lations , rotations  and  a change  of  scale.)  A general  account  of 
such  methods  is  given  by  Ackermann  (1972).  Although  technically 
classifiable  as  semi-analytical  methods  (by  virtue  of  the  instru- 
mentation formation  of  individiual  models),  they  differ  from  the 
semi-analytical  methods  of  Class  1 in  a significant  way,  namely, 
in  the  simultaneity  of  the  recovery  of  all  interconnecting  trans- 
formations. Because  of  this,  the  methods  of  Class  2 have  proven 
to  be  decidedly  more  effective  than  those  of  Class  1 and  accordingly 
merit  a separate  classification. 

In  the  United  States  the  expression  "simultaneous 


TABLE  2.7-1. 


CLASSIFICATION  OF  METHODS  OF  ANALYTICAL  AEROTRI ANGULATION 


CLASS 


1 


2 


3 


METHODS 


.Strip  Assembly  bv 
(Connective  Polv- 
' nomial  Transfor- 
mations 

Simultaneous  Adjust1 
'ment  of  Independent 
Models  by  Similar-  ! 
| ity  Transformat- 
ions 

Bundle  Adjustment 


[Bundle  Adjustment 
(with  Sel f-Cal ibra-  j 
l tion 


EXAMPLES  OF  COMPUTER  PROGRAMS 


NRC  System  (National  Research  Council,  j 
Canada)  USC&GS  Three  Photo  Relative 
Orientation 


ANBLOCK  (Vand  den  Haut) 

PAT-M  Series  (Ackermann,  et.  al., 

Stuttgart  University) 


BSOR-BLOCK  (DBA/RAOC/ACIC) 

USC&GS  Block  Aerotr iangulation  (NOS) 
MUSAT  (ETL/TOPOCOM) 

SURBAT  (DBA/RADC/ACIC) 

COMBAT  (DBA) 

I 

STEREO  (DBA) 

SMAC  (DRA/ETL) 

COMBAT  II  (DBA) 

BAP  (Bauer  & Muller,  Hannover ,W . G . ) 
Generalized  MUSAT (Automatic) 

STEREO  with  Geometric  Constraints (DBA) ' 


* This  listing  is  not  intended  to  be  all-inclusive. 
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block  aerotrinagulation : is  commonly  used  to  designate  the 
methods  of  Class  3;  Europeans  tend  to  favor  the  more  concise 
expression  "bundle  adjustment"  that  is  used  in  Table  1.  Methods 
of  Class  3 exercise  the  utmost  degree  of  mathematical  and  statis- 
tical rigor.  This  involves  the  simultaneous  least  squares  triangu- 
lation of  all  bundles  of  rays  from  all  exposure  stations  to  all 
measured  points  coupled  with  the  simultaneous  accomplishment  of 
the  resection  of  all  participating  photos  and  adjustment  of  the 
ground  survey. 

The  methods  of  Class  3 are  based  on  the  common  assumption 
that  measured  film  coordinates  have  been  adequately  corrected  for 
all  significant  sources  of  systematic  error  and  thus  are  contaminated 
only  by  random  error.  Those  of  Class  4 recoqnize  that  this  is  not 
necessarily  the  case  in  practice;  unresolved  systematic  errors  in 
film  coordinates  often  assume  appreciable  significance  and  must  be 
accorded  due  consideration  if  optimum  results  are  to  be  produced. 

This  can  be  accomplished  in  considerable  measure  by  a process  where- 
in additional  parameters  consisting  of  coefficients  of  appropriate 
error  models  are  recovered  within  the  bundle  adjustment  itself. 

Such  recovery  constitutes  a process  of  self-calibration  and  hence 
the  designation  "bundle  adjustment  with  self-calibration".  This 
last  catagory  of  programs  is  of  particular  importance  with  reqard 
to  improving  the  accuracy  of  point  positioning.  Two  programs  of 
specific  interest  are  the  Generalized  MUSAT  (Greve  et.al  1 976)  and 
the  program  STEREO  with  Geometric  constraints  (DBA,  1972).  These 
programs  contain  provisions  for  constraining  object  space  points 
to  lie  on  specific  geometric  surfaces  such  as  planes,  circular 


cylinders,  etc. 


In  addition  to  these  general  cateaor ios , special 
programs  have  been  developed  for  lunar  orbiter  and  Apollo 
Photography  which  incorporate  orbital  constraints  to  force  all 
exposure  stations  to  lie  along  an  arc  defined  by  six  continuous 
orbital  elements.  These  orbital  oarameters  replace  the  three 
exposure  station  position  parameters  associated  with  each  photo- 
graph and  add  six  parameters  per  strip.  These  additional  para- 
meters are  placed  in  the  border  of  the  normal  equation  (See 
Section  4.6  for  a detailed  development  of  the  banded-bordered 
normal  equations). 

DBA  Systems  has  also  developed  two  programs  for 
performing  analytical  trianculations  with  panoramic  photoqraphv. 
The  first  program  for  Rome  Air  Development  Center  (RADC)  included 
error  models  for  the  KC-56A  rotating  prism  camera  (RADC-tr-7 1-240) 
and  the  second  for  the  AdoI  lo  Optical  Bar  Panoramic  Camera 
(DBA,  1970).  These  programs  were  experimental  in  nature  and 
not  used  in  conventional  trianoulation . 


2 . 8 Coordinate  Determination 

As  indicated  in  the  previous  section  the  imaqe  space 
coordinates  are  a function  of  three  basic  classes  of  parameters: 
Exterior  orientation.  Systematic  deformation  parameters  and 
Object  Point.  Typically  the  interior  orientation  parameters  are 
calibrated  prior  to  a block  adjustment  and  are  considered  as 
known.  The  primary  purpose  of  the  block  adjustment  proqram  is 
to  derive  the  exterior  orientation  parameters  and  those  object 
points  carried  in  the  block  adjustment. 

There  are  frequently  many  points  not  carried  in  the 
block  adjustment  for  which  object  coordinates  are  desired.  These 
coordinates  are  determined  based  on  the  intersection  of  ravs 
from  two  or  more  photographic  imaqes  of  the  same  point.  The 
same  basic  colinearitv  equations  are  used  for  this  process  except 
that  the  adjustment  is  only  in  terms  of  the  object  coordinates 

(wv- 

It  should  also  be  noted  that  the  derived  positions 
are  in  a rectangular  coordinate  system,  therefore,  during  the 
last  phase  of  the  point  positioning  svstem  these  coordinates 
are  converted  to  the  desired  system  (UTM,  LSR,  PSR,  Geographic, 
etc . ) . 

2 . 9 Data  Base  Concept 

The  concept  of  a photoarammetr i c point  positioning 
data  base  is  illustrated  in  Figure  2.0-1,  assuming  the  material 
availability  question  is  answered  with  a ves  (In  Sections  2.1 
throuqh  2.8  it  was  assumed  the  answer  was  no). If  the  material  is 
available  and  if  it  has  been  previously  adjusted  it  is  only 
necessary  to  measure  points  on  the  photographs  and  perform  the 


-20- 


coordinate  determination  to  derive  object  space  position. 

2.10  Summary 

The  preceding  discussion  of  the  various  aerc- trian- 
gulation processes  are  not  meant  to  be  all  inclusive,  but  only 
to  serve  as  an  introduction  to  the  ensuing  analysis.  Section  3.0 
(Shortfall  Analysis)  provides  a brief  discussion  of  those  areas 
of  each  process  which  DBA  believes  may  need  improvement.  Section 
4.0  (Discussion  of  Alternatives)  describes  a number  of  alternative 
techniques  which  will  potentially  overcome  the  shortfalls  identified 
in  Section  3.0.  Lastly,  Section  5.0  (Conclusions  and  Recommendations) 
discussed  conclusions  derived  from  the  preceding  sections  and 


suggests  a systematic  approach  for  their  implementation. 


3.0 


SHORTFALL  ANALYSIS 


The  previous  section  outlined  the  current  state-of- 
the-art  in  analytical  aerotrianqulation . This  section  nresents 
a critical  review  of  this  system  and  identities  areas  where 
improvements  can  be  made  in  terms  of  (a)  Accuracy,  (b)  Time, 

(c)  Manual  Labor  involved  and  (d)  Skill  Level  required.  Each 
of  the  items  discussed  in  the  previous  section  will  be  analyzed 
for  potential  improvement  in  light  of  currentlv  available 
technologies . 

3 . 1 Material  Availability 

Current  available  material  determinations  are  directly 
related  to  an  individuals  familiarity  with  past  activities  in  an 
area.  The  amount  of  material,  as  well  as  quality  may  be  directly 
related  to  the  past  activities  of  only  one  organization  in  an 
area  and  may  not  consider  material  available  from  other  sources. 
Requirements  may  dictate  the  amount  of  time  available  for  material 
research  and  consequently  limit  the  magnitude  of  the  material 
search . 

The  basic  shortfall  in  determining  material  availability 
is  the  current  inability  to  present  to  any  user  at  any  time  all  of 
the  data  pertaining  to  a specific  area,  in  a timely  manner. 

3 . 2 Coverage  Planning 

Coveraqe  planning, as  currently  employed,  is  generally 
a manual  process.  As  with  most  manual  operations,  it  is  a tedious 
and  time  consuming  effort.  Consideration  must  be  given  to  such 
parameters  as  scale  of  the  final  oroject,  C-factor  of  the  eomoilation 
instrument,  type  of  acquisition  sensor  and  platform,  size  and  shape 


of  the  area  of  interest,  etc.  The  time  consumed  during  this 
process  drastically  reduces  response  tine  of  the  total  Aero- 
triangulation  System. 

3 . 3 Material  Acquisition 

The  process  of  material  acquisition  is  directly 
related  to  the  type  of  sensor  used,  incorporated  with  the  type 
of  platform  available  and,  as  previously  mentioned,  the  processing 
of  such  material.  Again  the  process  is  time  consuninq  and  related 
to  the  technologies  of  equipment  and  support  materials  rather 
than  being  an  analytical  problem.  As  such,  DB.n  has  no  intent 
to  delve  further  into  this  process. 

3.4  Image  Point  Selection , Identification  and  Mpnsurat ion 
Durinq  this  process,  a target  is  identified  through 

photographic  interpretation  techniques,  this  step  is  normal  1'’ 
performed  manually  by  specially  trained  photointerprotors . 
Additionally,  a target  mav  be  selected  and  identified  on  one 
type  of  photography  ar d must  be  transferred  very  precisely  onto 
a photographic  image  format  from  which  measurements  can  be  made. 
The  point  transfer  and  mensuration  steps  are  normally  performed 
manually  with  the  aid  of  electro-optical  measuring  equipment. 

The  total  process  is  extremely  time  consunina  as  a result  of 
the  manual  efforts  required. 
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There  are  three  basic  areas  where  steps  can  be  taken 


in  the  data  preprocessing  which  will  minimize  the  computer  run 
time  for  subsequent  block  adjustments.  These  areas  are: 

1.  Computation  of  initial  aporoximations , 

2.  Data  editing,  and 

3.  Data  arrangement. 

Obviously  if  initial  values  of  the  parameters  to 
be  adjusted  are  close  to  the  true  values  fewer  iterations  of  the 
least  squares  process  wil  be  required  (hence  less  run  time) . 

The  deviation  of  initial  approximations  for  some  parameters 
(such  as  exposure  station  positions)  can  be  very  tedious  and  is 
often  a guess  at  best.  To  rectify  this  problem,  techniques  for 
computationallv  deriving  initial  approximations  related  to  special 
types  of  photography  have  been  implemented  on  computer  assisted 
stereocomparators.  However,  these  techniques  have  not  been  applies 
to  aerial  photography  nor  to  monoscopic  instrumentation  or  as  a 
step  in  the  preprocessing  for  manual  comparators. 

The  second  area  ? data  editinq,  has  been  only  super- 
ficially addressed  for  aerial  photography . As  with  determination 
of  initial  values,  techniques  for  ore-adjustment  data  editinq  have 
been  implemented  on  a computer  assisted  storeocomputer . These 
techniques  could  be  modified  and  implemented  on  computer  assisted 
monoscopic  comparators  or  in  a preprocessing  program. 

The  last  area  which  is  insufficiently  exnlored  in 
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the  data  preprocessing  step  is  the  rearrangement  of  the  data  for 
efficient  formation  and  solution  of  the  normal  equations.  This 
step  in  most  cases  is  currently  performed  manually.  This  step 
is  slow  and  tedious  and  frequently  overlooked,  resulting  in 
excessively  long  run  times  for  jobs  which  could  be  performed 
much  more  efficiently.  This  is  particularly  true  with  non- 
parallel flight  lines  and  photography  of  widely  varying  scale. 

In  analyzing  the  automatic  reordering  techniques 
currently  employed  (Glaser,  1971),  it  is  clear  that  they  each 
suffer  from  one  or  more  of  the  following  limitations: 


1.  they  are  not  completely  universal  in  that  they 
reouire  a certain  degree  of  regularity  in  the 
block  and  cannot  efficiently  accommodate  photos 
of  greatly  varying  scales  or  in  randomly  oriented 
cross-strios ; 

2.  they  are  not  completely  automatic  in  that  they 
require  human  intervention  to  perform  a manual 
reordering  of  the  input  data  or  to  assess  the 
improvement  or  optimization  achieved; 

3.  they  employ  algorithms  that  do  not  provide,  either 
manually  or  automatically,  an  optimum  minimization 

4.  they  are  desianed  onlv  for  the  more  common 
systems  of  banded  normal  equations  and  cannot 
accommodate  the  more  complex  banded-  bordered 
systems . 


3 . 6 Block  Adjustments 

Current  analytical  adjustment  techniques  require 
improvements  to  facilitate  computations  in  three  broad  areas. 

The  first  area  is  the  excessive  computer  run  time  generated  with 
current  algorithms  due  to  inefficient  data  handling  algorithms 
or  improper  exploitation  of  existina  algorithms.  Another 
major  contributing  factor  is  the  incomplete  data  preprocessing 
(Sect  ion  3.5). 
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The  second  major  shortfall  of  current  aerotriangulation 
block  adjustment  is  that  current  formulations  provide  for  limited 
accuracy.  They  do  not  adequately  model  systematic  errors  nor 
do  they  consider  all  the  auxiliary  information  that  is  available. 

The  third  major  drawback  is  the  high  skill  level 
required  to  perform  a block  adjustment.  This  is  because  of  the 
complexity  of  the  algorithms  employed. 

Each  of  the  specific  programs  described  in  Section  2.7 
have  limitations  on  the  number  of  photographs,  number  and  tvpe  of 
error  models  constraints  etc.  This  class  of  limitation  is 
dependent  upon  the  specific  aloorithms  implemented  in  the  computa- 
tional equations.  Specific  techniques  used  by  individual  programs 
will  be  discussed  in  Section  4.0  as  general  methods  which  could 
be  implemented  in  other  programs  to  increase  their  capability. 

3 . 7 Coordinate  Determination 

The  speed  and  accuracy  with  which  point  coordinates 
may  be  determined  is  directly  related  to  the  speed  and  accuracy 
of  the  Block  Adjustment  and  Triangulation  Process.  The  limiting 
factor  for  this  process  is  the  turnaround  time  for  Data  Reduction 
at  a specific  facility.  Most  often, such  factors  as  priorities 
and  management  enter  into  this  process.  The  discussion  of  such 
factors  is  beyond  the  scope  of  this  report. 

3 . 8 Summary 

The  preceding  paragraphs  identified  specific  areas 

in  the  Analytical  Aerotriangulation  System  which  renuirc 
investigation  for  improvement.  Shortfalls  have  been  identified 
resulting  from:  inefficient  data  handling,  time  consuming  reduction 


models,  and  required  manual  processes.  Some  of  these  shortfalls 
have  been  identified  as  being  beyond  the  scope  of  this  report  or 
currently  under  investigation  through  other  RADC  efforts  and  will 
not  be  discussed  further.  However,  several  problem  areas  have 
been  identified  as  being  worthy  of  discussion;  Material  Availability, 
Coverage  Planning,  Data  Preprocessing,  and  Block  Adjustment.  Detailed 
discussions  for  each  of  these  areas  relating  to  alternative  techniques 
are  addressed  in  Section  4 of  this  reoort. 
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4.0 


DISCUSSION  OF  ALTERNATIVES 


Section  3.0  identifies  those  processes  in  the  Aero- 
triangulation  System  which  can  be  improved  through  the  use  of 
state-of-the-art  computation  technologv  and  advanced  algorithms. 
The  specific  areas  of  interest  include:  Material  Availability, 
Coverage  Planning,  Auxiliary  Sensors,  Ground  Control  Data,  Data 
Preprocessing,  and  Block  Adjustment.  Alternative  approaches 
and  techniques  for  each  of  these  are  presented  in  the  following 
subsect  ions . 

4 . 1 Material  Availability 

The  use  of  modern  mini-computers  combined  with  state- 
of-the-art  data  base  technology  can  be  applied  to  the  problem 
of  accessing  the  available  material  pertinent  to  an  area.  This 
approach  would  be  an  operator  interactive  system  which  would 
allow  the  user  to  interogate  a data  base  containing  information 
relating  to  the  coverage  available  over  an  area.  This  data  base 
would  include  algorithms  which  specify  the  status  and  quality  of 
this  coverage.  The  system  would  consist  of  a series  of  CRT 
stations  communicating  with  a master  data  base  facility. 

Such  a system  has  been  implemented  by  the  USGS  and 
is  currently  operational.  This  system  provides  the  user  with  a 
listing  which  identifies  the  tvoe  of  imacrery,  duality,  percent  of 
cloud  cover  data,  first  frame  center,  scale,  altitude,  and  limits 
of  coverage.  This  system  is  located  at  the  EROS  Data  Center  in 
Sioux  Falls,  South  Dakota.  This  master  data  center  is  accessed 
by  remote  terminals  at  several  USGS  sites. 

This  system  as  currentlv  configured  will  not  meet 


the  requirements  for  RADC  or  DMA,  however,  it  is  a base  for 
development  of  a system  tailored  to  specific  Air  Force  require- 
ments. 

4 . 2 Coverage  Planning 

The  manual  process  of  Coverage  Planning  can  be  improved 
by  use  of  a graphics  display  interactive  terminal.  Such  a system 
could  display  a map  of  the  area  to  be  photographed.  The  operator 
could  then  indicate  to  the  computer  the  approximate  location  of 
control  along  with  the  characteristics  of  the  sensor.  The 
computer  would  then  display  a flight  profile  including  a ground 
footprint  for  the  planned  mission.  The  system  would  then  perform 
a simulated  block  adjustment  and  display  the  theoretical  error 
contours  to  be  generated  by  such  a fliaht  plan.  The  user  would 
then  adjust  his  coverage  parameters,  control  specification  or 
auxiliary  data  requirements  to  optimize  the  coveraqe  so  that 
maximum  utilization  of  the  photogranhv  could  be  obtained  with 
minimal  expense. 

The  envisioned  system  could  then  produce  a hard  copy 
of  the  flight  plan  for  use  by  the  aircraft  navigator  in  acquiring 
the  photographs. 

Another  use  of  this  system  would  be  to  assist  strike 
mission  planners  in  developing  flight  paths  between  two  points 
by  including  such  thinqs  as  terrain  avoidance  algorithms.  For 
these  advanced  applications,  a digital  data  base  of  the  area 
would  be  required. 


4 . 3 


h. 


Auxiliary  Sensors  Onboard  Aircraft 
There  aro  ; nunber  of  auxiliary  sensors  which  are 
either  currently  available  or  will  be  available  in  the  near 
future.  These  sensors  have  the  capabilty  of  providing  an  external 
source  of  aircraft  position  and/or  attitude  which  could  be  used 
as  a supplement  to  the  ohotooraphic  system.  These  auxiliary 
sensors  provide  two  potential  benefits  for  aerotrianaulation . 

The  first  is  an  increase  in  accuracy  and  the  second  is  a set  of 
consistent  initial  approximations  to  the  exterior  orientation 
parameters.  The  ability  to  have  a consistent  set  of  initial 
approximations  is  valuable  in  that  preadjustment  editina  of 
data  is  more  reliable  and  the  number  of  iterations  of  the  least 
squares  algorithms  could  also  be  reduced. 

It  should  be  noted  that  many  of  these  auxiliary 
sensors  provide  consistent  but  systematically  biased  data.  These 
systematic  bases  can, in  some  cases  be  modeled  so  that  useful  infor- 
mation can  be  gained  from  the  sensors. 

The  Statoscope  and  Airborne  Profile  Recorder  are 
examples  of  such  auxiliary  sensors .The  indications  of  absolute  altitude 
and  absolute  change  in  altitude  are  of  insufficient  accuracy  to 
improve  the  accuracy  of  the  overall  trianciulation  due  to  system- 
atic biases  in  the  system.  In  general,  the  difference  between 
photoqranmetrieal ly  derived  exposure  stations  and  Statoscope 
altitude  observations  have  a systematic  set  of  residuals  which 
can  easily  be  modeled  in  the  block  adjustment.  These  parameters 
would  be  included  in  the  border  of  the  normal  equations  as  they 
would  be  common  to  a large  set  of  photographs  (Section  4.6.3 
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presents  a detailed  discussion  of  the  border  md  Section  4. 6. 1.4 
outlines  the  observational  ecruations  qenerated  by  the  use  of 
Statoscope  type  data). 

The  second  tyoe  of  sensor  is  the  hiqh  dynamic  Global 
Positioning  System  (GPS).  Tins  system,  as  currently  configured, 
uses  simultaneous  ranoe  observation  of  tour  satellites  to  compute 
positional  data  at  any  instant  in  time.  There  are  twenty  four 
satellites  in  the  proposed  GPS  conf 1 aurat ion  of  which  from  five 
to  twelve  will  be  visible  at  any  time.  The  current  desian  goal 
for  positioning  with  the  GPS  high  dynamics  system  are (S.\M80, 1974) 

Horizontal  8 Meters  CFP 

Vertical  10  vcters  PF 

Velocity  0.1  Knots  PE 

There  are  three  alternative  techniques  available  for 
use  of  this  data.  The  first  is  the  brute  force  approach  of 
inputting  the  observed  aircraft  position  at  the  time  of  the 
exposure  along  with  the  accuracy  estimates  indicated  above.  In 
most  cases  of  interest,  data  with  this  level  accuracy  will  add 
very  little  to  the  overall  accuracy  of  the  block  adjustment. 
However,  as  initial  approximations,  this  data  is  about  three  to 
four  times  more  accurate  than  that  currently  used  for  standard 
triangulation.  This  aspect  alone  would  bo  beneficial  from  the 
preadjustment  data  editing  and  possible  reduction  of  the  number 
of  iterations  required.  This  technique  will  require  no  change 
to  existing  trianoul ation  urograms  other  than  possibly  adding  a 
coordinate  transformation  step  in  the  data  processing. 

The  second  method  for  using  GPS  positional  data  is 

a 
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through  error  modeling.  The  errors  identified  above  will  not  be 


independent  between  exposure  stations.  That  is,  there  will  be 
systematic  biases  between  the  GPS  derived  positions  and  the 
photogrammetrically  derived  positions.  A large  portion  of  this 
error  can  be  attributed  to  the  ephemeris  errors  and  the  geometry 
of  the  specific  satellite  constellation  beina  observed  bv  the 
GPS  receiver  unit.  This  gives  rise  to  the  error  model  approach 
for  using  the  GPS  derived  position  as  a oart  of  the  photogrammetric 
reduction.  This  error  model  technique  emplovs  the  development 
of  a mathematical  model  v/hich  describes  the  difference  between 
the  GPS  positions  and  the  photogramraetricallv  derived  position. 

This  error  model  would  contain  a few  parameters  which  would  be 
common  to  a large  subset  of  photographs.  These  parameters 
would  then  be  added  to  the  border  of  the  overall  least  squares 
adjustment.  A mathematical  model  for  this  is  discussed  in 
Section  4.6.1. 3. 

A third  alternative  technique  actually  models  the 
continuous  motion  of  the  aircraft  in  soace  using  a spline  technique. 
The  GPS  receiver  can  be  modified  to  output  positional  data  at 
a very  high  frequency.  For  purposes  of  this  discussion,  assume 
that  GPS  is  capable  of  providing  data  at  a rate  of  thirty  times 
per  second.  This  data  would  be  used  to  determine  the  aircraft's 
position  in  space  with  the  photogrammetric  data  beinq  constrained 
to  lie  along  this  derived  function  wilh  appropriate  error 
model  terms  tc  model  the  systematic  discrepancy  between  the 
photogrammetric  and  GPS  positions.  The  actual  model  for  this 
system  is  discussed  in  Section  4. 6. 2. 5. 


J 
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Another  source  of  external  data  is  an  onboard  attitude 


and  attitude  rate  sensor.  This  data  can  be  used  nuch  as  was  done 
with  the  GPS  data  discussed  above.  Since  this  data  can  be  used 
by  the  sane  basic  techniques  as  the  GPS  data  (i.e.  direct  input, 
error  modeling  or  functional  models),  no  attempt  has  been  made  in 
this  report  to  explicitly  discuss  external  attitude  data.  However, 
one  should  bear  in  mind  that  the  techniques  discussed  for  the  GPS 
data  can  be  applied  to  externally  derived  attitude  data. 

4 . 4 Ground  Control  Data 

The  normal  type  of  Ground  Control  Data  used  for  block 
adjustment  is  the  position  of  photoidentif iable  objects.  This 
control  data  can  be  in  a variety  of  coordinate  systems  and  derived 
from  a number  of  sources.  Conventional  ground  surveying  is  the 
most  common  source.  However,  in  recent  years,  satellite  doppler 
surveying  has  become  an  extremely  attractive  method  for  determining 
precise  ground  positions  particularly  when  one  considers  that  this 
system  derives  points  relative  to  the  reference  ellipsoid, not  the 
Geoid  as  is  the  case  in  normal  surveying.  This  permits  one  to 
easily  develop  a function  relatina  the  local  Geoid  to  the  reference 
ellipsoid  in  the  area  of  interest.  This  is  accomplished  by 
observing  a few  of  the  conventional  surveyed  points  with  doppler 
tracking  devices  such  as  JMR  or  GEOCIEVER  and  by  empirically 
modeling  the  differences  between  the  derived  coordinates. 

In  addition  to  the  three  dimensional  absolute  control, 
there  are  other  types  of  information  which  could  he  exploited  in 
a block  adjustment  to  increase  the  accuracy  without  expensive 
surveying.  These  are:  (1)  relative  positioning  constraints, 
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(2)  inertial  survey  data  and  (1)  qeometric  constraints.  Each  of 
these  items  are  discussed  briefly  in  the  followinq  sections. 

4.4. 1 Relative  Pos i t i onina  Constraints 

There  are  a number  of  sources  of  control  information 
which  have  been  classified  as  relative  oositioning  constraints, 
four  of  which  are  listed  below: 

1.  Local  Survov 

2.  Distance  Data 

3.  Azimuth  or  Angle  Data 

4.  Equal  or  Differential  Elevation  Data 

Local  survey  data  consists  of  point  coordinates  which 
have  been  determined  very  precisely  bv  local  survey  but  have  not 
been  tied  into  an  absolute  datum.  Many  such  local  survevs  mav 
be  available  within  the  area  covered  by  the  adjustment.  These 
surveys  are  locally  consistent, but  are  biased  relative  to  each 
other  and  to  a master  datum.  These  biases  can  be  accommodated 
by  the  incorporation  of  datum  shift  parameters  in  the  border  of 
the  block  adjustment.  PLODS  (Haag  & Hodge  1972)  is  the  only 
program  to  our  knowledge  which  currently  has  provisions  for 
accommodating  these  datum  shifts. 

The  second  type  of  control  information  which  is 
available  is  the  precisely  measured  distance  between  two  points, 
the  locations  of  which  mav  be  unknown.  This  type  of  data  could 
be  collected  by  using  Geodiometer  type  distance  measurements 
across  larqe  lakes  or  between  mountain  tons  in  remote  areas. 

PLODS  and  SAPGO  have  provisions  to  accommodate  distance  constraints. 
SAPGO  requires  that  both  points  be  on  a sinqle  photograph.  This 
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restriction  is  not  in  PLODS  since  the  distance  constraint  is 
carried  in  the  border.  Section  4. 6. 2.1  describes  the  mathematics 
for  the  rigorous  inclusion  of  distance  constraints  in  the  border 
of  a conventional  block  adjustment. 

The  azimuth  of  a line  is  another  source  of  relative 
control  which  has  heretofore  not  been  properly  exploited  in 
analytical  aerotriangulations . As  with  distance  constraints 
both  PLODS  and  SAPGO  have  the  ability  to  handle  azimuth  constraints. 
SAPGO  is  limited  in  the  single  Photo  constraint  whereas  PLODS 
utilizes  the  border  for  azimuth  constraints.  The  mathematics 
for  adding  the  azimuth  constraints  to  an  aerotriangulat ion 
reduction  is  presented  in  Section  4. 6. 2. 2 

The  equal  and/or  differential  elevation  constraints 
are  also  available  in  SAPGO  and  PLODS.  DBA  has  also  included 
this  capability  in  our  proprietory  block  ad justment-COMBAT-II . 

This  constraint  is  generated  by  the,  fact,  that  frequently  points 
are  known  to  have  equal  elevations,  however,  the  value  of  that 
elevation  is  unknown.  An  example  of  this  type  of  information 
is  the  points  along  the  shoreline  of  a lake.  The  method  by  which 
this  data  has  been  included  in  PLODS  and  COMBAT- I I is  described 
in  Section  4. 6. 2. 3 of  this  report. 

4.4.2  Inertial  Survey  Data 

Within  the  past  two  years,  new  tvnes  of  qround  survey 
system  have  appeared  on  the  market. These  are  the  inertial  surveying 
systems  such  as  the  Position  and  Azimuth  Determination  Systems  (PADS) 
and  the  Litton  Autosurvevor . These  svstems  utilize  acceleration 
observations  derived  from  three  orthogonal  qyroscopes.  These 
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accelerations  are  nummericallv  integrated  with  respect  to  tine 
by  an  onboard  computer  to  produce  velocity  and  to  provide 
position  data.  This  integration  is  performed  every  seventeen 
milliseconds . 

The  normal  operation  of  the  inertial  surveying  unit 
is  to  derive  standard  point  positioning  information.  This  type 
of  object  point  data  can  be  directly  entered  into  an  analytical 
triangulation  as  a control  point.  The  accuracy  of  this  data  would 
be  marginal  in  that  the  following  are  the  svstem  specifications 
for  the  (PADS) (ETL, 1975); 


Horizontal 

20 

Me te  r s 

CEP 

Vertical 

10 

Meters 

PE 

Azimuth 

1 

Mil 

RMS 

However,  accuracies  of  1 meter  horizontal  CEP  and  .5  meters  Vertical 
(PE)  have  been  demonstrated  with  the  Litton  Autosurveyor  by  Gregarson 
in  CANADA  (Gregarson,  1975).  This,  however,  is  after  adjustments 
have  been  made  to  known  control  with  a systematic  error  distribution 
function  applied.  The  basic  model  used  by  Gregarson  is  of  the 
following  form: 

AF  - A<j>  = AAA  + BA  At  + EA$  + 0 

AL  - AA  = AA<{>  + BA<t>t  + EAA  + D 

where  F,L  are  the  true  latitude  and  longitude 

4>,A  are  the  measured  latitude  and  longitude 

A-  is  the  alignment  error 

B-  is  the  drift  in  alignment 

E-  scale  error  of  quantizer 

D-  error  accumulator 

t-  time  relative  to  start  of  tranverse 
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The  Gregarson  parameters  (A,  B,  E,  D)  could  easily  be  included  in 
the  analytical  model  for  aerotrianoulation. 

Another  approach  to  including  inertial  surveying 
observations  directly  in  an  aerotrianoulation  is  through  the  use 
of  an  autoregressive  procedure.  This  technique  has  been  included 
in  DBA's  Doppler  Survey  reduction  program  and  was  documented  by 
Brown  and  Trotter  in  1969.  However,  the  use  of  the  autoregressive 
model  in  aerotriangulation  has  not  been  fully  exploited.  The 
general  concept  for  its  implementation  is  presented  in  Section  4.6.4. 

Another  attractive  method  for  exploiting  the  capability 
of  the  inertial  survey  system  is  in  conjunction  with  the  geometric 
constraints  discussed  in  the  next  section. 

4.4.3  Geo met ric  Constraints 

In  the  past,  analytical  aerotriangulation  has  used 
only  specific  object  points  as  the  basic  reference.  That  is,  the 
same  objective  point  is  measured  on  multiple  photographs.  In  recent 
years  programs  have  been  developed  which  incorporate  some  specific 
relationships  between  points  such  as  distance,  azimuth  and  equal 
elevation  constraints  as  discussed  in  Section  4.4.1. 

Greve  et  al  in  1976  describes  a program  called  the 
"GeneraLizt  I Multi~Phot  Analytical  Bl  ■'<  Adjustment  Program”-  This  program 
was  designed  to  handle  a large  number  of  non-metric  type  photographic 
views  of  the  same  general  object.  This  program  includes  algorithms 
for  constrainting  specified  sets  of  object  points  to  lie  on  one  or 
more  geometric  surfaces.  The  types  of  surfaces  which  this  proqram 
can  accommodate  are: 

• co-circular  points 
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co-suher ical  ooints 


• co-cylindrical  points 

• co-linear  points 

• points  lying  on  parallel  lines 

• coplaner  ooints 

• points  on  two  planes 

• points  at  equal  intervals  on  a line 
Properly  exploited  geometric  constraints  may  be  of 

tremendous  practical  utility  in  the  analytical  reduction  of 
digital  photographs.  This  can  be  demonstrated  by  a simple  hypo- 
thetical example.  Assume  we  have  two  aerial  photographs  of  a road. 
This  road  is  known  to  be  straight.  Then  a]l  points  measured  along 
this  road  should  be  constrained  to  line  on  a straight  line  in  object 
space.  Assume  further  that  the  edae  of  the  road  is  digitized  on 
each  photograph  independently.  This  implies  that  the  points  measured 
on  one  photograph  have  no  direct  correspondence  to  ooints  measured 
on  the  other  photoqraph.  In  other  words  all  points  measured  on 
each  photograph  must  be  considered  from  a classical  viewooint  as 

single  rav  points.  In  classical  aerotriangulation  the  only  single 

* 

ray  points  which  contribute  to  the  solution  are  control  points. 

It  is  highly  desirable  to  develon  algorithms  wherebv  these  two 
sets  of  single  ray  points  could  be  constrained  to  lie  along  a straight 
line  in  object  space.  A proposed  method  for  developing  such  an 
algorithm  is  discussed  in  Section  4. 6. 2. 4. 

In  general,  there  are  a number  of  geometric  constraints 
of  this  type  which  can  be  exploited  in  analytical  reductions.  A 

* Single  ray  pass  points  are  freouentlv  carried  in  the  solutions 
when  sub-blocks  are  reduced  for  editing  purposes.  However,  those 
points  are  weighted  so  that  they  do  not  affect  the  solution. 
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typical  constraints  arc: 


• Sections  of  Roads  and  railroads  can  be  defined  as; 

Horizontal  straioht  lines  or, 

* 

Horizontal  sections  of  a circular  curve. 

Vertical  straight  lines  or, 

Vertical  sections  of  parabolic  curve. 

• Intersecting  roads  can  be  constrained  to  the 
same  elevation. 

• Points  lyinq  on  the  perimeter  of  objects  which 
have  known  ueometric  shapes  (such  as  storaqe 
towers)  can  be  constrained  to  lie  on  the  appropriate 
mathemat ical  surface. 

The  use  of  this  type  of  constraint  will  become  more  important  as 
laroer  focal  length  imaaing  systems  come  into  use,  particularly  with 
regard  to  the  digital  imaging  sensors. 

4 . 5 Pat  a Preprocess i no 

As  stated  in  Section  2.6,  Data  Preprocessing  can  be 
divided  into  four  basic  functions:  Transformation,  Generation, 

Editing  and  Arrangement.  Section  3.6  identified  the  shortfalls 
associated  with  the  last  three  areas.  In  this  section  alternative 
toheniques  are  out  1 ined,  wh ich»  if  implemented  could  possibly  alleviate 
some  of  these  shortfalls. 

■1 .1.1  Generation  of  Initial  Approx imat  ions 

Tn  classical  bundle  adjustment  aerot rianaulat ion 
processes,  there  are  two  classes  of  parameters  which  must  bo 
approximated.  These  are  the  elements  of  exterior  orientation 
and  the  pound  point  coordinates  for  all  pass  points  (Tnit  i il 

* Spiral  curves  are  sometimes  used  for  turnpike  int<  i chanqos 
but  we  can  omit  them  from  practical  cons i dot  at i on  at  t Ip 
present 
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values  for  control  points  are  set  equal  to  their  observed  value.) 


As  mentioned  earlier,  the  problem  of  determining 
initial  approximations  to  the  exterior  orientation  is  of  primary 
importance.  An  attractive  approach  to  this  can  be  envisioned  if 
one  realizes  that  the  normal  equation  is  rather  insensitive  to 
the  actual  measurements  and  are  nuite  sensitive  to  the  general 


location  on  the  photograph  and  the  number  of  points  on  the  photo- 
graph. The  classical  normal  equation  can  be  represented  as 

6 = N-1  C = (BTWB)  "1B'rWr  4. 5. 1-1 

In  general,  both  the  weight  matrix  (w)  and  the  iacobian  matrix 

(B)  are  insensitive  to  the  actual  observations,  but  the  discrepancy 

vector  is  very  sensitive  to  the  observation.  It  one  realizes  this, 

T -IT 

the  (B  WB)  B W portion  of  eauation  4. 5.1-1  could  be  precomputed 

based  on  nominal  values  for  the  parameters  and  a nominal  image 

point  distribution.  Multiplying  the  initial  discrepancy  (•)  bv 
T -IT 

the  matrix  (B  WB)  B W would  then  result  in  a set  of  corrections 
to  the  nominal  values.  These  corrections  would  then  be  applied  to 
the  initial  approximations  for  the  riaorous  least  squares  solution. 
It  should  be  noted  that  if  one  considers  the  weight  matrix  as  a 
constant  diagonal  matrix  with  enual  diagonal  elements,  it  can 
be  factored  out.  The  net  result  is  that  this  nominal  matrix  is 
exactly  the  same  as  the  "Desian  matrix"  employed  by  Rauhala (1976)  . 
Rauhala  has  demonstrated  that  this  desian  matrix,  when  inter- 
preted properly,  can  be  used  to  simplify  some  least  squares  problems 
and  yield  the  same  results  as  normal  techniaues.  The  use  of  this 
method  is  very  attractive  for  computing  initial  exposure  station 
values . 
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A practical  method  of  implementing  this  technique  could 
be  the  generation  of  a theoretical  design  matrix  using  9 points  in 
classical  locations.  As  the  points  are  read  into  the  preprocessor, 
the  9 points  closest  to  those  used  in  the  generation  of  the  theor- 
etical design  matrix  could  be  used  to  compute  the  initial  exterior 
orientation  approximations.  Expanding  this  technique  to  derive 
object  space  coordinates  may  also  be  possible. 

Another  alternative  is  the  use  of  a recursive  algorithm 
for  the  computation  of  initial  approximations.  DBA  has  implemented 
such  a system, on  the  TP-3/P1  computer  assisted  stereocomputer  for 
a special  type  of  photography  . 

In  general  terms  this  sane  algorithm  would  be  implemented 
in  an  off-line  data  preprocessing  program  or  on  a monoscopic  computer 
assisted  comparator.  The  basic  philosophy  of  this  recursive  tech- 
nique is  that, each  time  a point  is  read  (either  on  the  comparator 
or  into  a computer),  its  contribution  to  the  exterior  orientation  is 
computed.  This  contribution  is  then  tested  (data  editing)  prior 
to  adding  it  to  the  previous  estimate  of  the  initial  value.  There- 
fore, as  each  image  point  is  input,  the  inital  values  are  refined. 

The  result  of  this  step  is  that  the  initial  parameter  values  which 
are  input  to  the  block  adjustment  are  much  closer  to  the  true  values 
and  the  data  is  more  unlikely  to  contain  blunders. 

4.5.2  Data  Editing 

Data  Editing  can  be  divided  into  two  basic  categories  - 
Manual  and  Automatic.  Manual  being  that  the  user  phvsically  detects 
and  eliminates  the  aberrant  observation.  Automatic  means  that  the 
computational  algorithms  detect  and  handle  these  blunders. 
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However,  all  automatic  editinq  must  be  subject  to  manual  review. 

4 . 5 . 2 . 1 .Manual  Editing  Algorithms 

The  efficiency  and  to  some  extent  the  ability  of  an 
operator  to  manually  edit  data  is  directlv  related  to  the  manner 
and  amount  of  data  presented  to  him.  An  attractive  approach  for 
manual  editing  involves  the  use  of  a CRT  display  interacting 
with  a computer  which  is  capable  of  editing  the  data  and  nictor- 
ially  presenting  the  results  of  the  editing. 

There  are  basically  two  techniques  for  implementing 
this  type  of  system.  The  first  is  through  the  use  of  a remote 
terminal  attached  to  a large  scale  CPC  which  is  also  used  to 
perform  the  least  squares  adjustment.  The  second  technique  is 
through  the  use  of  a mini-computer  as  a dedicated  editing 
station. 

In  either  case,  the  results  of  the  preprocessing  could 
be  tabulated  and/or  graphically  displayed.  These  displays  should 
be  optimized  for  convenience  of  the  user  from  an  analysis  stand- 
point. That  is,  all  of  the  results  should  be  available  for 
display.  However,  key  items  should  be  displayed  so  that  the  user 
is  not  required  to  go  through  many  pages  of  numeric  tabulations 
to  detect  a few  blunders.  This  implies  that  the  routines  should 
have  some  type  of  blunder  detection  algorithms  so  that  potential 
blunders  can  be  displayed  in  a manner  which  is  readily  apparent 
to  the  user. 

4 . 5 . 2 . 2 Automatic  Edi  ting 

During  the  preprocessing  step  the  major  problem  with 
automatic  editing  of  pass  points  is  that  image  measurement  residuals 
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are  a function  of  the  moan  object  space  coordinates  which  in  turn 
are  a function  of  the  image  measurements.  This  inter-relationship 
makes  is  extremely  difficult  to  detect  blunders.  A solution  to  this 
is  available  when  points  appear  on  note  than  two  photographs.  That 
is,  if  a point  appears  on  three  or  more  photoqraphs  in  a block 
there  are  two  or  more, two  ray  combinations  of  intersections  which 
can  be  analyzed.  A pair-wise  analysis  of  the  intersections  of  the 
rays  will  provide  editing  insight  into  the  duality  of  the  other 
rays  for  this  point.  This  analysis  can  be  in  terms  of  residuals 
or  the  computed  object  space  coordinates. 

Another  potential  technique  for  preprocessing  editina 
is  an  analysis  of  the  difference  between  observed  object  space 
coordinates  of  control  points  and  their  computed  value  based  on 
the  initial  exterior  orientation  approximations.  These  differences 
could  be  used  as  the  basis  for  determining  systematic  discrepancies 
in  the  approximations  to  the  exterior  orientation  parameters. 

Yet  another  potential  technique  for  editing  data  and 
computing  the  initial  approximations  to  the  exterior  orientation 
parameters  is  a point  by  point  updating  scheme,  whereby  each  point 
is  input  to  the  program  and  its  contribution  to  the  exterior  orien- 
tation parameters  is  computed.  This  change  in  orientation  is  then 
used  as  a criteria  for  editina  (See  Section  4 . 5 . i ) . 

4.5.3  Reordering 

Section  2.6.1,  discussed  current  techniques  for 
reordering  the  parameters  to  bo  adjusted  in  an  analytical  aero- 
tr iangulation  to  minimize  the  bandwidth  of  the  normal  equations 
thereby  increasina  the  speed  and  efficiency  of  the  solution.  The 


basic  short  comings  of  these  techniques  were  then  summarized 
in  Section  3.6.2.  In  sumary,  current  techniques  are  limited 
for  one  or  more  of  the  following  reasons: 

1)  they  are  not  completely  universal, 

2)  they  are  not  completely  automatic, 

3)  they  do  not  orovide  ootimum  minimization,  or 

4)  they  can  not  accommodate  banded  bordered  systems. 

The  most  advanced  of  the  current  photogrammetric  re- 
ordering or  optimization  techniques  is  that  used  in  program  COMBAT. 
This  advanced  algorithm  was  discovered  and  adapted  to  photogrammetric 
application  as  part  of  an  in-house  investigation  performed  by  DBA 
in  1971.  At  that  time,  a study  was  made  of  the  latest  research  in 
bandwidth  minimization  being  conducted  in  other  fields,  such  as 
structural  stress  analysis,  which  also  encounter  scarce  diaaonally 
banded  matrices.  Several  promisina  techniaues  (Elphingstone , 1 972 ) , 
(Akyuz  & Utku,  1968),  (Rosen,  1968)  were  adapted  to  the  ohoto- 
arammetric  application  and  evaluated  for  speed  and  correctness 
of  minimization.  The  most  promising  of  these  was  (Glazer  & Saliba, 
1971)  subsequently  incorporated  into  Proqram  COMBAT. 

Although  significant  improvements  are  yet  possible  in 
the  speed  and  efficiency  of  this  algorith,  its  primary  limitation 
is  its  current  inability  to  handle  banded-bordered  systems. 

As  regards  future  development  and  improvement  of 
photogrammetr ic  bandwidth  minimization,  it  is  highly  recommended 
that  continuing  research  of  this  be  sponsored, 
objectives  of  this  research  would  be  to: 
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The  primary 


4 . 6 


Block  Adjustments 


k 


A number  of  techniques  are  available  for  increasina 
the  efficiency  and  accuracy  of  photoqrammetr ic  block  adjustments. 

These  have  been  grouped  under  the  following  broad  headinqs  for 
purposes  of  this  discussion: 

1)  Error  Modelinq 

2)  Constraints 

3)  Banded-Bordered  Structure 

4)  Autoregressive  Modelinq 

5)  Automatic  Editing 

The  first  two  items  include  methods  whereby  additional 
information  can  be  added  to  the  standard  photogrammetric  reduction. 

Many  of  these  methods  result  in  a set  of  normal  equations  which 
have  a " banded-bordered"  structure.  An  efficient  alqorithm  for 
solving  special  types  of  banded-bordered  structures  is  present'-' 
in  Section  4.6.3. 

Section  4.6.4  describes  an  approach  whereby  serially 
correlated  observations  can  be  rigorously  incorporated  in  a least 
squares  adjustment.  As  will  bo  shown, this  technique  provides 
an  alternative  method  for  handling  some  error  sources  which  have 
heretofore  been  modeled  by  empirically  derived  polynomial  or 

Fourier  series.  j 

The  last  area  whore  major  technological  emphasis 
is  required  is  in  the  automatic  editinq  aspects  of  a photogrammetric 
block  adjustment.  Each  of  these  Five  major  areas  are  described 
in  detail  in  the  following  paragraphs. 
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4.6.1  Error  Mode  1 i ng 

Systematic  errors  have  been  classically  handled  in 
photoqrammetrv  b''  applvinn  corrections  to  the  basic  observations 
in  the  oreprocess ina  nhase.  However,  for  ultra-precise  photo- 
grammetric  tr iangulation  and  for  using  observational  data  qathered 
bv  some  tvoes  of  auxiliarv  sensors  it  is  necessarv  to  carry 
parameters  which  model  systematic  errors  as  an  integral  nart 
of  the  triangulation  reduction.  Section  4. 6. 1.1  reviews  the 
general  concept  of  including  systematic  error  model  parameters 
in  a photoqrammetr ic  block  adjustment.  'this  is  followed  in 
Section  4.6. 1.2  throuoh  4.6. 1.4  with  the  development  of  specific 
error  models  for  Anomalous  Distortions,  Ops  Positioning  errors 
and  Statoscope  Observational  errors. 

4-6.  1.1  Model inq Systematic  errors  in  a photoqrammetr  ic 

reduction . 

4 . 6 . 1 . 1 . 1 General  Formulation 

The  classical  photoqrammetr ic  bundle  adjustment  is 
based  upon  the  colinearity  condition  equations  which  can  be 
represented  functional ly*as 


where 


r 1 T 


: fc(Xi' 


X . 
1 


C) 


4 . 6 . 1 - 1 


i - represents  information  pertaining 

to  the  i-th  photograph 

represents  information  pertaining 
to  the  j-th  object  point. 


' 

| x 


denote  the  image  space  coordinates  of 
of  the  j-th  point  on  the  i-th  photograph. 


* Explicit  definition  of  this  set  of  equations  can  be  found  in 
the  Manual  of  Photoqrammetrv. 
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xc 

1 

- 

is  a three  element  vector  of  exposure 
station  position  for  the  i -th  photograph. 

1 <f>c 

1 

- 

is  a three  element  vector  of  orientation 
angles  for  the  i -th  photograph. 

1 j 

- 

is  a three  element  vector  of  object  space 
coordinates  for  the  j-th  around  point: 

c 

- 

is  a three  element  vector  of  interior 
orientation  parameters 

The  vector  C in  the  above  representation  contains  the  camera 
focal  length,  and  the  x,y  principal  point  coordinates.  These 
elements  are  normally  treated  as  constants  in  a photogrammetric 
reduction  , hence  they  are  separated  from  the  parameters  by  a semi- 
colon . 

The  above  expression  is  valid  ,only  if  the  true  value 
of  all  the  parameters  and  the  image  coordinates  are  known,  and 
that  all  systematic  errors  have  been  properly  accounted  for. 


In  reality  observations  for  the  image  coordinates 

are  subject  to  random  error.  By  denoting  observations 

'o' 

by  a superscript  this  relationship  can  be  expressed  as: 

[x,y]T.  = [x  , y]  °T  + [v  , v ] T . 4 . 6 . 1 - 2 

L L ' * J ij  Lx  yJii 

where  the  v's  denote  observational  residuals.  Substituting  this 
expression  into  4. 6. 1-1  results  in  the  following  model. 
oT  . r iT 


[x,y]  . . + [v  ,v  ] . • = f ( X . , f.,  X . ; C) 
L ' J l ] L x y J l j l ' l ' i 


4.6. 1-3 


This  represents  the  classical  observational  equations  used 

c c 

in  bundle  adjustments  where  the  vector  X , <,*>.  and  X.  represent 

l i t 

the  parameters  to  be  determined  by  a least  squares  adjustment. 
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However,  if  there  are  uncompensated  systematic  errors, 
equation  4. 6. 1-2  is  invalid  and  should  be  rewritten  as  follows. 

[x.yli'j  = + tvx  vy]T  + g (p)  4. 6.1-4 

where  g denotes  a systematic  error  model  with  parameters  P. 

The  exact  form  of  the  error  model  depends  on  the  physics  of  the 
system.  (Specific  examples  will  be  demonstrated  in  Section  4. 6. 1.2 
through  4. 6. 1.4.  ) The  parameter  vector  P will  contain  the  p 
parameters  required  to  model  the  error.  Substituting  this  new 
representation  of  the  true  coordinates  into  the  orojective 
equation  results  in 

[X,y]  i j +K  Vylij  + 9(P)  =f(Xi'  VXj;C>  4. 6. 1-5 


which  can  be  rewritten  as 
i OT 


[*,y]°j  +K  vy]ii  = h(xi-  V xj'  P:C) 


4. 6. 1-6 


Equation  4. 6. 1-6  now  represents  the  basic  photogrammetric 

observational  equations. 

Assuming  that  h is  a non-linear  function  of  all 
c c 

parameters  (X^,  $ ^ , X^  , P) , an  iterative  technique  is  required 

to  solve  for  the  parameters.  We  shall  linearize  the  basic 

observational  equations  about  some  arbitrary  initial  values  of 

the  parameters  and  about  the  observed  image  measurements.  The 

true  value  for  the  parameters  are  qiven  in  equation  4.6. 1-7,  where 
* oo  * 

the  superscript  denotes  the  arbitrary  initial  approximations 

and  <5  represents  corrections  to  these  aonrox  imations . 
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Xc  = [xc]°°  + *x 


iJ 

r.cioo 


4. 6. 1-7  a 


, c _ r , c I uu  ip 

= l*  J + 6 1 


X.  =[X.]  + 6X 

3 L 3J  3 


P = [P]°°  + 6 


Substituting  equations  4. 6. 1-7  into  4. 6. 1-6  results  in 


r 1 OT  r 
I x , y I . . + I 


oo 


V V 


c oo 


>■  x y-*  1 3 " i X ^ 


h (X  . +fvc  , 1 +4  7 , XT  + ,P  + ;c) 


.‘1  l x • 


4. 6. 1-8 

or  by  simply  rearranging  terms  the  following  function  can  be 


developed . 


oo 


0=F(x°.+v  , y.  . +v  ,XC  + 6 vc  , 

l]  x . . 1 l -j  y..  3 X. 

j 13  j r13  j 3 


:C°  +•  c ,X°°  +6  ,POCV  ;c) 
3 4 ■ 3 x p 

3 3 


4 . 6 . 1 -9 


Equation  4.6.1-  9 can  nQw  be  approximated  by  the  first  two 
terms  of  the  Taylor  series.  This  linearization  process  results 
in  the  following  matrix  equation. 


• c c . c c 

A . . V . f BX  . A x +B : . v!  + P . P + B . - 3 x . + I . . =0  4.6.1-10 


13  13  13  3 13  3 13 


13  13  13 


where  the  following  matrix  definitions  hold. 

A . 

(2,2) lj 


pF.  . / 

I i3/'x..,y 


/ i 3 i 3 
T 


4.6.1-11  a 


(2,1) 


13 


[v  v ] 

x yJ 
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4.6.1-11  c 


x 

B.  . 
1J 
(2,3) 

, c 


B 

(2,3) 


-3 


8.  . 

13 

( 2 , p) 


B.  . 
13 

(2,3) 

(2,1) 


ij/3X^  ] 

ij/H^j 

ij/3P 


-[■ 


ij/3X. 


o o 

X . . , V . . 

1 3 13 


oo 


F (X: 


y°°  p°°. 

' 3 ' ' 


C) 


d 


e 


f 


<7 


f*  , >7*  , 6^,  arc  defined  by  equation  4. 6. 1-7. 

The  dimensions  under  the  matrices  are  defined  by  the  fact 
that  Fisa  two  element  vector  function  (one  Function  for  x and 
one  for  y)  which  results  in  the  A and  B matrices  having  2 rows. 
The  number  of  columns  in  each  Jacobian  matrix  is  determined  by 
the  number  of  elements  in  the  parameter  vector,  i.e.  X*? , ( 

X_.  each  contain  3 elements  whereas  P contains  p elements. 

Each  of  these  partial  derivative  matrices  is  evaluated  at  the 
initial  approximation  to  the  parameters.  The  discrepency  vector 
t is  the  difference  between  the  observed  image  coordinates  and 
the  function  F evaluated  at  the  approximations. 

Note  that  equation  4.6.1-10  is  linear  with  respect 
to  the  parameter  changes  ( d ) . Equation  4.6.1-10  represents  the 
linearized  observation  equations  for  the  j-th  object  point  imaged 
on  the  i-th  photograph.  Assume  that  there  are  a total  of  m 
photographs  and  a total  of  n points  in  the  photoqrammet ric 


block  to  be  adjusted. 
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Also  assume  for  the  moment  that  all  points  are  imaged  on  all 

■k 

photographs  . First  we  shall  collect  all  the  linearized 
observational  equations  generated  by  the  j-th  ground  point. 


This  set  of  equations  is 


_>r  x 
♦ 3ij  °i 


.5 . 1-12 


* 3b 


The  offset  between  the  A's,  B ' s,etc.  is  easily  seen  if  one 
considers  that  the  partial  derivatives  of  the  image  points  on 
the  i-th  photograph  with  respect  to  the  parameters  for  the( i + I) st 
photograph  are  zero.  In  general  the  derivatives  are  zero  for 
all  photographs  other  than  the  i-th  photograph. 

Equation  4.6.1-12  can  be  rewritten  more  compactly 
using  matrix  notations  as 

A.  V.  + B*C  6xC  + B*C  5*C  + §,«P  +B.<tX  + . = 0 

13  1 1 1111 

4.6.1-13 


where  the  following  definitions  are  employed 


A.  . 

i ll 


4.6.1-14  a 


A . = 
1 

( 2m, 2m) 


A . I 

m]  J 


This  second  assumption  simplifies  the  presentation  of  the 
mathematical  development  and  will  be  dropped  later  in  thi: 
r eport . 


_ to 


(2m, 3) 


4.6.1-14  h 


lj 

(2m, 1) 


p x 

6 and  6 ^ are  defined  in  equation  4.6.1-7d  and  c respectively. 

The  next  step  is  the  collection  of  all  linearized 


observation  equations  for  the  n object  points.  This  set  of  equations 


As  with  equation  4.6.1-12,  this  set  of  equations  can  be 
rewritten  in  matrix  form  as  : 


E = 

( 2mn , 1 ) 


e 2 
£ 3 


4.6.1-17  h 


6X  and  a1*1  are  defined  by  equations  4.6.1-14  e and  f 

respectively.  6^  is  defined  by  equation  4.6.1-7d.  The  elements 

of  the  other  matrices  in  equation  4.6.1-17  are  all  defined  by 

equation  4.6.1-14.  Note  that  there  are  2mn  equations  and 

c c 

2mn  + 6 m + 3n  + p unknowns  (v's,  6X  's,  6'  's,  6x's,  and  SP,s) . 

Obviously  there  are  more  unknowns  than  there  are  equations  so 
that  there  are  an  infinite  number  of  solutions  which  will  satisify 
equation  4.6.1-16.  The  most  desirable  solution  is  the  one  which 
satisfies  the  linearized  observation  equation  and  simultaneously 
minimizes  the  sum  of  the  squares  of  the  weiqhted  residuals,  i.e. 

T 

s - v w v — > minimum  value  4.6.1-18 


It  has  been  demonstrated  by  Brown  1955  that  this  minimization  is 
accomplished  by  solving  the  following  equations: 


N 6 = c 


where 


T T - 1 

B (AAA  ) B 

T T -1 

B (AAA  ) . 


4.6.1-19 


4.6. 1-20 


4.6. 1-2] 
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4.6.1-22 


4.6. 1-23 


The  covariance  matrix  (A)  and  weight  matrix  (w)  are  discussed 
in  depth  in  Section  4. 6. 1.1. 4 and  the  solution  of  equation  4.6.1-19 
for  the  parameters  6 is  discussed  in  Section  4.6.3. 

Before  proceeding  with  the  specific  examoles  for  error 
models,  there  are  three  additional  areas  which  must  be  discussed. 
These  are;  (1)  the  rearrangement  of  parameters , (2)  the  introduction 
of  parameter  observations  and  (3)  the  weighting  of  observations. 
These  areas  are  described  in  the  following  subsections. 

4. 6. 1.1. 2 Rearrangement  of  Parameters 

In  eauations  4.6.1-11  through  4.6.1-23  the  exoosure 

. c 

station  (X  ) position  was  carried  seoaratelv  from  the  orientation 
angle  parameters  ($  ).  Since  these  two  classes  of  oarameters  are 
associated  with  a single  photogranh,  it  is  freguently  convenient 
to  define  a six  element  vector  of  exnosure  station  oarameters  (E) 
which  is 


4.6.1-24 
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Substituting  4.6.1-24  into  4.6.1-10,  with  the  usual  definition 
of  initial  approximations,  results  in  the  following  set  of  linearized 
observational  equations  for  the  i-th  point  imaged  on  the  i-th 

photograph 

. . oo  .... 

A . . v . . + B..6.  + B . . o + B . . 6 . . + e . .=  0 4.6.1-25 

l]  13  13  1 13  13  13  13 

where 


and  all  other  elements  are  identical  to  those  defined  bv  eauation 
4.6.1-11.  The  collection  process  proceeds  as  was  described 
for  equations  4.6.1-12  through  4.6.1-17  with  the  obvious  definition 
substitutions.  The  net  result  is  a linearized  set  of  observational 

equations  which  can  be  represented  in  matrix  terms  as: 

. . 00 

Av  + B<5  + BS  + B6  + e = 0 4.6.1-26 

This  is  the  general  form  of  linearized  observational  equations 
which  one  would  normally  utilize  to  represent  a photogrammetric 
triangulation  with  systematic  error  modeling.  However, in  other 
cases, it  is  desirable  to  separate  the  position  elements  from  the 
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orientation  elements  as  exoressed  by  eouation  4.6.1-16.  During 
the  discussions  which  follow  in  Sections  4. 6. 1.2  through  4.6. 1.4 
specific  examples,  will  be  disscussed,  some  of  which  can  be 
better  treated  using  the  form  represented  by  4.6.1-16  and  others 
by  the  form  4.6.1-26. 
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4.6.1. 1.3  Parameter  Observations 


As  discussed  in  Section  4. 6. 1.1.1  the  observational 
equations  are  linearized  about  arbitrary  initial  approximations. 
It  is  reasonable  to  assume  that  in  some  cases  these  parameters 
can  be  directly  or  indirectly  observed.  The  purpose  of  this 
section  is  to  outline  a procedure  for  including  observational 
equations  for  parameters  directly  into  the  overall  system  of 
linearized  observational  equations  prior  to  initiating  the 
least  squares  adjustment. 


The  inclusion  of  directly  observed  parameters  is 
a straight  forward  process  which  has  been  documented  in  many 
reports  (Ref.  Brown, et  al  1965 , Davis&Riding : 1 970 , Haag&Hodqe  1972) 
In  general, it  proceeds  as  follows.  Assume  that  observations 
for  the  ground  point  coordinates  are  available.  The  true  value 


for  these  coordinates  is 


where: 


X . 
1 


X°  + vx 
1 1 


4 . 6 .1-27 


is  the  three  element  vector  renresentino  the 
true  position  of  the  object  ooint. 

X°  denotes  the  observed  value  for  these  object 

1 'c 

space  coordinates.  (Notice  the  superscript 

is  used  to  denote  observed  quantities) 
x 

Vb  is  the  ihree  element  vector  of  residuals. 

Equation  4.6.1-27  can  be  substituted  directly  into  equation 
4.6.1-7c  as  follows- 


X°  + vx 

1 1 


x°°  + 


which  by  appropriately  rearranging  terms  becomes 


,.x  Xx  , x _ _ 

V.  - 6 . + e =0 

1 1 


4.6.1-28 


where 


ex  = x°  - x°° 

1 1 


There  will  be  one  such  equation  generated  for 
each  observed  object  space  point,  it  is  assumed  herein  that 
all  n object  points  are  observed(this  assumption  will  be  discarded 
later).  The  collection  of  all  sets  of  equation  4.6.1-28  results 
in  the  following  matrix  equation; 


„x  X . X _ _ 

V -8  + e =0 


4.6.1-29 


where 


VX  = 
( 3n , 1 ) 


V 


| V 

L n 


X 

e = 


( 3n , 1 ) 


X 

n 


X 

and  8 is  identical  to  the  definition  in  4.6.1-17q. 

The  logical  question  arises  as  to  why  not  simply 

define  the  arbitrary  initial  parameter  approximation  (x”°)  as 

the  observed  value?  This  obviously  reduces  the  iX  vector  to 
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zero.  However,  one  must  remember  that  the  parameters  will  be 
iteratively  solved.  At  the  completion  of  each  solution  the 
6 vector  is  applied  to  the  previous  parameter  approximations  and 
the  solution  is  again  performed.  This  correction  cannot  be  applied 
to  the  observation,  therefore,  after  the  first  iteration  the 
approximation  used  for  the  linerization  processes  will  be  different 

x 

from  the  observation,  hence  the  t vector  will  be  non-zero. 

Equation  4.6.1-27  implies  that  only  random  errors  are 
associated  with  the  observation  of  the  j-th  point.  This  may  not 
be  the  case.  Just  as  with  image  measurement  there  may  be  systematic 
errors  associated  with  the  observation  of  the  j-th  point.  These 
systematic  errors  will  be  common  to  many  other  (if  not  all)  object 
space  points.  For  each  discussion  it  is  assumed  that  the 
systematic  errors  are  in  fact  common  to  all  points  and  can  be 
modeled  by  a function  q'  which  has  a set  of  parameters  P'  just 
as  for  the  imaqe  point  observations. 

Equation  4.6.1-27  then  becomes: 

X . = X°  + VX  + q'  (P’  ) 4.6.1-30 

1 1 1 

which  when  substituted  into  equation  4.6.1-7c  becomes: 

X°  + VX  + q ' ( P'  ) = X°°  +<SX  4.6.1-31 

1 3 3 3 
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Proceeding  as  before  this  results  in  a linearized 
set  of  observational  equations  for  the  j-th  points  which  are 


X 0 > P 
V + B.  6 

3 1 


6 + e . = 0 

3 3 


4.6.1-32 


This  set  of  equations  can  be  directly  added  to  the  structure 
given  by  4.6.1-16  as  outlined  below.  First  redefine  the  vector 
P to  include  the  parameter  P'  i.e. 


P = 


4.6. 1-33 


The  partial  derivative  matrix  B. . must  also  be  redefined  to 

13 

include  the  partials  of  the  image  measurement  function  (F)  with 
respect  to  the  parameters  P1 . These  partials  will  be  zero.  So 
that  B becomes 


B..  J;,P^ 
13 


/ :'P 


•f  n r *f 
x3/  »p - 1 =!  ii/ap0 


/:,pj  -L 


4.6.1-34 


The  reverse  is  true  of  B^'  it  must  be  augmented  by  the  derivative 
of  the  object  point  equation  with  respect  to  the  image  point 
error  model  parameter  (P)  i.e.  ; 


o 

B’. 

3 


'9-s 


7 .p 


V p' 


3 


/3P 


I 


4.6. 1-35 


Using  this  new  definition  of  B all  object  point  observational 
equations  can  be  collected  so  that  the  following  matrix  representation 
is  valid. 


..  o o 

V+Bi'*-<5+f=0 


4.6.1-36 


L 


where 


r x i 

£i  i 


(3m, 1) 


( 3m, p+p 1 ) 


,X  X 

o , e 

m J [ m 


5 is  the  error  model  parameter  correction  vectors  defined  by; 


o 


n ?7_ 


4.6.1-37 


Rewriting  both  equation  4.6.1-16  and  4.6.1-36 
the  matrix  definitive  for  B in  equation  4.6.1-20  becomes: 

Xc  XC  <tc  ic  00  "" 

V + B 5 + Bv  6*  + B5+B6  + e = 0 

..  oo- 

V + B' 6-6+e  = 0 4.6.1-38 


o 

B B 

o 

B'  -I 


The  e vector  in  equation  4.6.1-21  must  be  augmented  by  t.  It 
should  be  noted  here  that  the  B matrix  has  been  redefined 
by  equation  4.6.1-34. 


Throughout  this  section  consideration  has  been 
given  to  observations  of  object  space  points.  This  is  presented 
only  as  an  example  of  the  technique  by  which  direct  observations 
of  parameters  can  be  included  in  a photogrammet ric  adjustment. 
In  fact  all  parameters  are  treated  just  as  the  object 
points  discussed  herein. The  net  result  is  a matrix  equation 
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representing  all  observations  as: 


where 


4.6.1-39 


In  this  case  the  error  model  parameter  vector  P is 


with  P'  being  the  error  model  parameters  associated  with 
the  image  points. 

xc 

P being  the  error  model  parameters  associated  with 
the  exposure  station  parameters. 

.c 

P’  being  the  error  model  parameters  associated  with 

the  orientation  angles,  nnd; 

X 

P being  the  error  model  parameters  associated  with 
the  object  space  coordinates. 

Kquation  4.6.1-39  represents  the  general  linearized 
form  of  til'  rv ■> ' "I'-’  .-.Mint  rnn. 
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1.6. 1.1. 4 Weighting  of  Observations 

Throughout  this  section  have  been  discussions  of 
observations  and  the  making  of  assumptions  concerning  the 
observability  of  all  parameters.  All  observations  must  include 
an  estimate  of  their  accuracy.  The  technique  for  utilizing 
these  apriori  accuracy  estimates  in  a photogrammetric  triangu- 
lation are  discussed  in  this  paragraph.  A natural  tall-out  of 
this  discussion  is  a convenient  technique  tor  dealinq  with  the 
assumption  made  earlier. 

First  i the  manner  by  which  apriori  accuracy  estimates 
for  image  measurements  are  utilized  in  rigorous  triangulation 
programs  will  be  described.  It  shall  be  assumed  throughout  this 
section  that  all  observations  are  independent  (i.e.  the  measure- 
ment of  two  image  points  are  not  correlated.)  Representing  the 
appropriate  covariance  matrix  for  the  measurements  of  the  i-th 
object  point  on  the  i-th  photogranh  by  A the  normal  covariance 
matrix  definition  is  obtained. 


A . 
il 


where 


4.6.1-40 


2 


o 


X 


2 

o 


y 


denote  the  x and  y variance  of  the  imane 
measurements . 


a 


xv 


a 

yx 


denote  the  covariance  between  the  x and 
y coordinates. 
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Normally  the  x and  v accuracy  estimates  are  equal  and  uncorrelated. 

a2  = o2  and  o =0.  However,  if  transformations  are  applied 
x y xy 

to  the  image  measurement  prior  to  the  adjustment,  then  the 
accuracy  estimates  must  also  be  transformed.  In  many  cases  this 
transformation  of  coordinates  causes  the  entire  matrix  (4.6.1-40) 
to  fill  in.  For  this  reason,  the  apriori  covariance  matrix 
will  be  considered  as  a full  2x2  matrix  . 

Equation  4.6.1-19  indicated  that  the  desired  solution 
for  6 using  the  weighted  least  squares  involved  the  formation  of 
normal  equations  (N)  by  the  formula: 

N = B T (A  A AT)-1B  4.6.1-41 

If  A is  an  identity  matrix  (which  it  normallv  is  in  photoqrammetr i r 

T —1 

adjustments) then  N, becomes  BAB.  It  is  convenient  to  define  the 
weight  matrix  W as  the  inverse  apriori  covariance  matrix 
i . e . . 

I 

w.  . = A -1  4.6.1-42 

13  13 

In  Section  4. 6. 1.1.1  it  was  assumed  that  all  points  appeared  on  all 
photographs.  This  assumption  v;as  necessary  to  simplify  the 
collection  of  the  observational  equations  using  standard  matrix 
notations.  Obviously  this  assumption  is  not  true.  If  a point 


does  not  appear  on  a photograph,  then  the  observation  of 

this  point  can  be  considered  as  havinq  an  infinite  variance.  This 
implies  that  the  weight  matrix  associated  with  this  point  is 

zero.  Theoretically  one  can  consider  all  the  non-existent 

* 

observational  equations  as  having  a zero  weight  matrix. 

Just  as  was  done  for  the  observational  equations  ,the 
weight  matrices  for  all  image  points  viewing  a given  point  are 
collected  into  a matrix  Vb  which  is  ; 


4.6.1-43 


and  then  for  all  n object  space  points  as  ; 


4.6.1-44 


Computational  algorithms  do  not  employ  this  technique,  they 
simply  use  only  the  actual  observational  equations. 
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In  both  of  the  above  cases  the  independence  of  observations  is 
implied  by  the  use  of  zero  covariance  terms  relating  observations 
of  one  image  point  to  the  next.  Section  4.6.4  shall  describe  an 
algorithm  for  handling  serially  correlated  data  which  violates 
this  basic  assumption  of  observational  independence. 

As  described  in  Section  4. 6. 1.1. 2,  the  image  measurements 
are  not  the  only  source  of  observational  data  to  be  considered  in 
a photogrammetric  adjustment.  Since  parameters  are  considered  as 
observable  quantities,  they  also  have  associated  with  them  an 
apriori  uncertainty.  These  apriori  accuracy  esitimates  are  handled 
exactly  the  same  as  those  described  above  for  the  image  points. 

4. 6. 1.2  Anomalous  Distortion  Error  Model 

It  has  been  demonstrated  by  Brown,  1968,  that  after 
compensating  for  calibrated  radial  and  decentering  distortions 
in  aerial  photographs,  there  remain  systematic  patterns  to  the 
image  measurement  residuals  after  a least  squares  adjustment. 

This  pattern  is  often  not  readily  apparent  when  looking  at  the 
residuals  on  a single  photograph,  but  it  does  become  very  apparent 
when  the  residuals  for  a large  number  of  photoqraphs  are  plotted 
on  a single  frame.  The  exact  source  of  these  systematic  errors 
has  not  been  identified,  however,  one  can  assume  that  the  primary 
contributors  are  platten  unflatness  and  uncompensated  film  defor- 
mation. An  approach  to  modelinq  this  systematic  error  is  throuoh 
the  use  of  qeneral  polynomials.  Equation  4. 6. 1-4  states  that  the 
true  value  of  an  image  coordinate  is  equal  to  the  observation  plus 
a random  residual  plus  an  error  model  function  g. 
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To  model  anomalous  distortion  we  may  define  g as  follows: 


1 


q = a +a. x + a„y  + a,xy  + 

^x  o 1 2 3 

gy  = bo  + bj_x  b2y  + b3xy  + 

4.6.1-44 

The  a's  and  b's  are  the  parameters.  The  number  of  parameters 
required  will  depend  on  the  order  of  the  polynomial  emoloyed. 

The  example  described  by  Brown  in  1968  indicated  that  a fifth 
order  polynomial  was  the  upper  limit. 

Instead  of  the  polynomial  approach  out  '.ned  above 
it  may  be  desirable  to  define  the  error  model  function  as  a 
general  two  dimensional  Fourier  series. 

As  described  by  Brown  1974,  a model  of  platten 
unflattness  could  also  be  included  in  the  general  adjustment  by 
defining  g as  follows; 

gx  = F(ci*2+c2xy+c3y2  + c4w3+C5x2y+c6xy'-+c7y3) 

V 2 2 3 2 2 3 

qy  = r^clx  +c2xy+c3y  +C4X  +C5X  y+c6xy  +c?  ) j 

4.6.1-45  J 

In  this  case  the  polynomial  (with  parameters  models 

the  unflatness  of  the  platten.  Again  a Fourier  series  could  be 
used  in  the  place  of  the  polynomial. 

4 . 6 . 1 . 3 G.P.S.  Positioning  Model 

As  described  in  Section  4.3  the  Global  Positioning 
System  instrumentation  can  be  mounted  on  aircraft  to  orovide 


J 


2 2 

a ,x  + acy  + . . . 
4 5J 

b4x2  + b5y2  + . . . 
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aircraft  positioning  at  specific  time  intervals.  The  current 
expectations  are  that  this  positioning  data  will  have  an  apriori 
uncertainty  of  8 meters  horizontally  and  10  meters  vertically. 

This  level  of  accuracy  will  not  appreciably  increase  the  accuracy 
of  photogrammetric  triangulations.  However,  it  will  reduce  the 
time  for  performing  an  adjustment, because  it  provides  better 
initial  approximations  and  has  many  editing  benefits. 

These  accuracy  estimates  are  relative  to  an  absolute 
reference  system  and  a majority  of  the  error  budget  is  contributed 
by  uncertainties  in  the  orbital  ephemeris  of  the  GPS  satellites. 
Since  the  same  satellite  configuration  will  be  tracked  by  the 
GPS  receiver  for  a period  of  time,  these  orbital  uncertainties 
will  be  common  to  a large  number  of  photographs.  For  purposes 
of  this  discussion,  assume  that  the  same  satellites  are  tracked  for 
all  photographs  in  the  block.  The  exposure  station  coordinates 
can  then  be  represented  by  the  following  matrix  equations: 

o 

Xi  = Xi  + Vi'  + g (P)  4.6.1-46 

o 

where  the  observation  is  the  GPS  observed  position  at  the 

time  of  exposing  the  i-th  photograph.  The  error  model  function 
g(P)  will  model  the  systematic  differences  between  the  GPS  position 
and  the  photogrammetrically  determined  position.  The  exact  form  of 
this  function  can  only  be  determined  after  experimental  data 
has  been  collected  , however,  it  is  reasonable  to  expect  that  a 
low  order  polynomial  for  each  component  will  be  adequate. 

-73- 


1 


The  assumption  of  using  the  same  satellite 
constellation  throughout  an  entire  photographic  mission  is 
probably  not  reasonable.  This  can  be  accommodated  in  a basic 
model  by  a defined  series  of  error  model  functions  (g),  one 
for  each  constellation  selected.  The  specific  function  to  be 
applied  for  any  given  photograph  will  be  determined  by  the  time 
interval  over  which  the  constellation  was  tracked  and  the  time 
of  the  exposure.  This  can  be  thought  of  as  reinitializing  the 
error  model  whenever  the  GPS  receiver  changes  constellations. 

4. 6. 1.4  Statoscope-  Observations 

The  Statoscope  is  a device  which  measures  differential 
changes  in  elevation  of-  an  aircraft  in  flight.  This  instrument 
can  be  keyed  to  the  aerial  camera  so  that  height  differences 
are  recorded  at  each  exposure  station.  The  absolute  accuracy 
of  these  differences  depends  to  a large  extent  on  the  barometric 
characteristics  of  the  atmosphere.  Assume  that  an  aircraft  is 
flying  at  a constant  elevation  but  that  the  barometric  pressure 
is  chanqina  along  the  flight  oath.  The  statoscope  would  record 
this  as  a change  in  aircraft  height.  In  realty,  the  output  of 
the  statoscope  is  a function  of  barometric  pressure  chanqes  and 
aircraft  altitude  chances.  It  is  reasonable  to  assume  that  the 
contribution  from  chanqinq  barometric  surfaces  can  be  modeled 
using  a low  order  polynomial  (probably  an  offset  term  and  a drift 
term).  This  can  be  included  into  the  overall  scheme  by  defining 
the  error  model  function  for  the  exnosure  station  appropriately. 
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4.6.2 


Constraints 


k 


During  the  Data  Acquisition,  Reduction  and  Block  Adjust- 
ment Phases,  additional  parameters  relevant  to  object  space  param- 
eters may  be  made  available  throuqh  auxiliary  data.  Such  data  may 
infer  specific  tolerances  to  which  we  may  wish  to  constrain  or 
enforce  relative  point  parameters.  This  section  defines  the  develop- 
ment of  algorithms  for  the  inclusion  of  such  constraints  into  the 
normal  equations.  This  will  result  in  a stronger  solution  while 
at  the  same  time  maintaining  the  block  diagonal  integrity  of  the 
required  portions  of  the  coefficient  matrix. 

4.6.2. 1 Distance  Constraint 

The  distance  constraint  allows  one  to  enforce  a given 
tolerance  on  the  chord  distance  between  two  surface  ooints.  As 
an  example  let  the  USR  coordinates  of  the  first  point  be  represented 
by  (X^ , Y^,  Z ^ ) and  the  second  ooint  by  (X.,,  Y.,,  Z0).  Thus  the 
distance  constraint  mav  be  represented  by: 

vdU  = 'Z  (xrx2)2  + (YrV2  +~<vyr  (4.6.2.i-u 

Let  us  now  define  the  partial  derivatives  of  this  constraint 
with  respect  to  the  USR  coordinates  of  the  surface  point;  such 
that ; 


We  may  now  derive  the  observational  equation  for  distance  constraints. 
First  of  all  let  us  define: 


\ " 3(XrVrZ1)  ' c2  3(F.2,Y2,22)  ' 


(/1 .6 .2.1-4) 


. ')  oo 

£ = V,  - V, 

di2  d12 


in  which  v°  are  observed  parameters  and  \u°  are  anproximated  parameters. 
di2  U12 


Thus  the  observational  equations  arc; 


v ' = u 6 + u 6 - e' 

d Cj  1 c2  2 


(4. 6. 2. 1-5) 


which  may  be  written 

v’  = u6-e'  (4.6.2. 1-6) 

Let  us  now  assume  a typical  set  of  normal  equations  of  the  form; 

<5  = N_1c  (4.6.2.1-7a) 

which  may  be  expressed  through  partitioning  by  ; 


= 

tv) 

-T 

N 

-1 

c. 

<s 

N 

N 

c 

(4.6.2. 1-  7a) 


N N - i N 

where  fl  my  be  of  the  form  o and  M my  be  of  the  form  = 

N ‘ N ' I 


lb  introduce  the  distance  constraint  let  us  solve  for  N and  C from 
equation  (4.6.2. 1-  7a)  in  the  form  of; 

N = BT  WB 
T 

C = B W 


w 


such  that; 


6 = (uT  w*  u)_1  (uTw'e')  (4. 6. 2. 1-8) 

in  which  w'  is  the  corresponding  constraint  weight  matrix  and  equation 
(4. 6. 2. 1-8)  is  the  equation  to  be  introduced  into  equation  (4.6.2.]- 7b)  , and 
results  in: 


— 

-1 

4 

N N 

c 

6 

NT  (uTw'u)  + N 

.T 

(u  we  ) + c 

- 

As  distance  constraints  may  be  applied  to  any  combination  of  points,  the 
constraint  matrix  may  have  elements  located  in  an”  position,  thus  destroying 
the  bloc/  diaaonality  of  the  N matrix.  In  order  to  maintain  the  dinconelity 
of  the  N matrix,  for  folding  and  solution,  these  constraints  must  be  introduced 
via  an  alternative  method.  Let  us  introduce  a”  additional  set  oc  urt-nr-Ti 
parameters  4'  such  that  equation  (4. 6. 2. 1-9)  becomes  . 
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= 

-T  " i .T 
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4’ 

= 

0 u j-w' 

-1 

(4.6.2.1-10) 


Equation  (4.6.2.1-10)  may  now  be  arranged  as  follows  in  order  to  reqain  the 
structure  of  equation  ( 4. 6. 2.1- 7b) : 
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Thus, the  constraint  effects  are  placed  into  the  border  of  the  handed- 
bordered  system  and  the  banded  integrity  of  the  N mtrix  is  retained. 

4. 6. 2. 2 Azimuth  Constraint 

Another  type  of  constraint  which  one  may  wish  to  apply  to 
ground  point  data  is  the  Azimuth  Constraint.  Such  a constraint  states 
that  the  angle  between  a line  from  the  first  point  (X^,Y^,Z^)  to  a 
reference  (eg.  True  North)  and  a line  frctr.  the  first  point  (X.  ,Y.,Z.) 
to  a second  point  is  held  to  a given  tolerance.  This  azimuth 

constraint  my  be  represented  by; 


(4.6.2 .2--} ) 


The  partial  derivatives  of  this  constraint  may  now  be  defined  with  resfoct 
to  the  USR  coordinates  of  the  surface  points  such  that; 
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(4. 6. 2. 2-2) 


in  which; 
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.6. 2. 2-3) 


The  observational  elements  (u  , u , • ')  in  equation  (4.6.2.1-11)  my  now 

C1  c2 

lx?  defined  with  respect  to  the  azimuth  constraint  mra  eters.  Following  the 
identical  development  represented  inequations  (4.6.3. 1 s)  through 
(4.6.2.1-11)  it  can  bo  shown  tint  the  azimuth  constraints  can  also  !»■  placed 
in  the?  torder  structure  so  as  not  to  dearado  the  N structure. 
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4.6.2. 3 Elevation  Constraint 

One  additional  constraint  which  may  be  applied  to 
ground  points  is  the  elevation  constraint  wherein  the  elevation 
difference  between  two  points  is  maintained  to  a desired  tolerance. 
Let  us  assume  a first  point  with  geographic  coordinates  (<f> , A ^ , 1^ ) 
and  USR  coordinates  (X^,Y^,Z^)  and  a second  point  with  geographic 
coordinates  ( 4>  2 , A 2 , h2 ) and  USR  coordinates  (X2,Y2,Z2).  The  equal 
elevation  constraint  between  these  two  points  may  be  represented 
by; 


(4. 6. 2. 3-1) 


Again  we  define  the  partial  derivatives  of  this  constraint  with 
respect  to  the  USR  coordinates  to  obtain; 


8vh 

12 = (cos<}>  cosA.  , cos^t  sinA.  , simtq  ) (4.6.3. 3-2) 

J(X1,Y1,Z1)  1 11  1 


(cos.}>2  cos.\2  , cos4>2  sinA^  , simf^) 


Again,  the  observational  elements  (u  , u , > ')  are  defined  in 

C1  c2 

equation  (4. 6. 2. 1-4)  with  respect  to  the  equal  elevation  constraint 

parameters.  Similarly  the  development  represented  in  equations 
(4.6.2. 1-5)  through  4.6.2.1-11)  is  followed  to  place  these  con- 
straints in  the  border  structure  and  again  maintain  the  integrity 
of  the  N matrix. 
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4. 6. 2. 4 


Geometric  Constraints 


As  identified  in  Section  4 . 4 . 3,  a valuable  source 
of  information  on  aerial  photographs  is  the  geometric  constraints 
imposed  by  the  imagery.  There  are  two  basic  approaches  to  this 
problem.  The  first  is  the  employment  of  identified  object  space 
points  and  constraining  them  to  lie  along  specific  mathematical 
surfaces.  This  is  the  approach  used  by  Greve  et.al.  1976  and 
Brown  et.al.  1971. 

The  second  approach  to  include  Geometric  constraints 
is  based  upon  a technique  developed  bv  Brown,  1967  called  The 
Method  of  Continuous  Tvaoes.  The  techniaue  was  developed  for 
precise  geodetic  positioning  using  ballistic  cameras  tracking 
a satellite. 

The  basic  concept  is  illustrated  in  Figure  ^.6.  2-1. 
Three  ballistic  cameras  (1,2,3)  are  viewing  a single  satellite 
pass  (c) . On  each  photograph  a number  of  points  are  measured 
along  the  trace.  The  same  point  is  not  necessarily  measured  on 
multiple  photographs.  All  of  the  points  measured  on  a photograph 
generate  a ruled  surface  illustrated  by  S^,S2,  and  in  the 
figure.  These  three  surfaces  must  intersect  along  the  actual 
trajectory  of  the  satellite  (c) . In  this  application, the 
orientation  of  each  photograph  was  known  by  stellar  reduction 
and  it  was  assumed  that  two  of  the  stations  were  at  known  locations 
with  the  only  unknowns  being  the  exposure  station  coordinates 
of  the  third  station. 
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Brown  demonstrated  that  a single  arc  was  theoretically  sufficient 
for  complete  determination  of  the  point,  however,  this  arc  was 

only  theoretical  and  not  practical.  However,  if  multiple 
satellite  passes  were  observed  from  a given  station,  then, the  three 
dimensional  coordinates  could  be  determined.  The  solution  developed 
additional  strength,  as  one  would  expect,  when  orbital  constraints 
were  applied  to  the  curve  C. 

The  idea  of  using  the  continuous  trace  concept  for 
modeling  geometric  constraints  in  aerial  photogrammetry  is  very 
attractive.  The  overall  technique  can  best  be  understood  by 
considering  a simple  case  - namely  a straight  line  in  object 
space.  Assume  that  this  line  is  imaged  on  two  photographs. 

Further  assume  that  a series  of  measurements  are  made  along  the 
line  on  both  photographs.  In  theory  the  line  should  be  straight 
on  both  photographs.  The  measurements  along  the  line  on  one 
photograph  and  the  exposure  station  define  a plane  in  space  which 
must  contain  the  object  space  line.  Since  two  such  planes  are 
generated,  one  for  each  photograph,  then  their  intersection 
will  be  the  line  in  object  space.  Assuming  that  no  other  infor- 
mation is  available  , the  intercepts  of  this  line  cannot  be 
determined  but  its  slopes  can  be.  In  fact,  the  intersecting 
planes  will  contribute  directly  to  the  determination  of  camera 
attitude,  but  not  position.  The  additional  fact  that  straight 
object  space  lines  should  be  straight  irrtaqe  space  lines  will  also 
strengthen  the  recovery  of  some  of  the  calibration  or  error  model 
terms  discussed  in  Section  4.6.1. 
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The  concept  of  intersecting  ruled  surfaces  can  be 
expanded  to  include  not  only  planes  but  generalized  surfaces. 

The  advantage  of  the  ruled  surface  concept  over 
those  heretofore  implemented,  lies  in  the  ability  to  measure  any 
number  of  points  along  an  image  on  an  individual  photograph  with- 
out transferring  specific  points  to  overlapnina  photographs. 

Further  advantages  of  this  type  of  constraint  may  be 
realized  through  the  use  of  slope  data  derived  from  an  inertial 
surveying  system  such  as  PADS. 

This  type  solution  would  also  make  better  use  of  the 
available  imagery  on  a photograph  than  is  done  by  the  classical 
point  triangulation  techniques. 

4. 6. 2. 5 Spline  Constraints 

The  use  of  spline  functions  in  photogrammetr ic 
triangulation  to  represent  error  model  functions  is  very 
desirable  in  that  a number  of  low  order  polynomials  can  be  used 
to  represent  a highly  erratic  function.  This  is  accomplished  by 
forcing  the  polynomials  to  be  continuous  at  the  end  point  and 
to  have  equal  first  derivatives  at  the  end  points. 

This  type  of  function  is  potentially  very  useful 
in  modeling  the  actual  movement  of  an  aircraft  as  a function  of 
time.  This  type  of  constraint  is  an  alternative  method  for 
incorporating  the  GPS  data.  Section  4.6.1  described  a method 
whereby  the  GPS  position  at  the  time  of  each  exposure  could  be 
included  by  modeling  the  difference  between  the 
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GPS  position  and  the  photoqrammetrically  derived  position  with  a 

polynomial.  In  this  section  a orocedure  shall  be  outlined  whereby 

the  actual  motion  of  the  aircraft  can  be  modeled. 

This  can  best  be  explained  b'f  considering  a single 

strip  of  photography.  The  time  for  acquiring  this  strip  will 

start  at  t and  continue  until  time  t,  . This  time  interval  will 
o k 

be  divided  into  k segments  with  end  points  and  center  points  as 
illustrated  below. 


?o  ^ h t2  *2  h *3  t4  *4  H t6  S 


FIGURE  4. 6. 2-2 


A series  of  third  order  oolynomials  will  be  used  to 
represent  the  position  of  the  aircraft  during  each  time  span 
(t{  1<t<t?M  =1,2, 3... k).  The  midpoint  of  each  time  span 
is  t , ( i = l , 2 . . . k) . For  each  time  span  there  will  be  three  polynomials 
(one  for  each  coordinate  direction) . For  the  X coordinate  the 
following  family  of  polynomials  will  be  used. 
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4. 6. 2. 5-1 


sf  (t)  - «£  a ♦ a^(t-tk)2  a tk_iit<lk 

In  order  to  insure  that  the  model  is  not  discontinuous  at  the 
end  points  we  shall  enforce  the  constraints  that  the  function 

S and  S must  have  the  same  value  at  the  common  end  point . Further , 

l £+1 

the  aircraft  velocity  must  be  the  same  at  that  end  point.  These 
two  constraints  can  be  specified  as: 

End  point  continuity 


S1  (tl}  = S2  (tl) 

s2  (t2)  = s3  (t2) 


4.6.2 .5-2 


Sk-1  (t  k-ir\  ^k-l1 


w 


Velocity  Continuity 
si  (V  = s2  (t  _1) 

s2  (t2)  = s3  (t2) 


4. 6. 2. 5-3 


Sk-1 (tk-l) 


Sk  (tk-l} 


where 

S , denotes  the  i -th  polynomial  in  any  of  the 
coordinate  directions  (Xc,  Yc,  Zc) 

S , denotes  the  time  derivative  of  the  f-th 

polynomial  in  any  of  the  coordinate  directions 


mhe  process  of  enforcing  equality  of  higher  order  time  derivatives 
could  also  included  if  required.  However,  for  our  purposes 

the  equal  velocity  constraint  at  the  end  points  should  be 
adequate.  In  general  terms  the  S-th  polynomial  can  be  written 
as 


)+ae3(t‘ 
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elsewhere . 
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4.6.2. 5-4 
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This  equation  is 
the  vector  Xc. 
our  colinearity 


l*,y 


']!: 


substituted  directly  into 
This  results  in  the  basic 
equations  being: 
fp(S(t.) , x ; c) 


equation 

functiona 


1 


4 . 6 . 1-  1 for 
form  of 


4. 6. 2. 5-5 


where  S represents  the  entire  family  of  polynomial  functions 


and  t^denotes  the  time  of  exposing  the  i-th  photograph. 

Xc  Xc  Xc 

Naturally  S is  a function  of  the  parameter  a ^ a 2 a 3 
xc  ? c 

a,  , a^, , . . . . a'  . The  direct  observations  of  aircraft 
<4  1 1 44 

position  from  the  GPS  receiver  are  entered  directly  by  the 
observation  equations  of  the  form  4.6.1-30  where  the  junction 

XC 

g(P)  is  replaced  with  the  family  of  polynomials  S where  B 

is  defined  as  the  derivative  of  S with  respect  to  the  parameters 

a^  (4  = 1,2,3  k;  s = 1,2, 3,4;  9 = X°  , YC , ZC)  . 

4s 

Two  sets  of  "dummy"  observations  are  introduced  to 
facilitate  the  incorporation  of  the  spline  continuity  constraints. 
These  dummy  observation  equations  are: 


p*  = se-l(Vl'  - Wl> 


D ? = si-i(VlJ  - Wl> 


e = 2,3 k 

4. 6. 2. 5-6 

We  are  not  really  observing  P or  D,but  because  we  desire  to 
enforce  the  continuity  constraints  thes_  values  should  be  zero. 

We  will  therefore  assign  observations  to  them  of  zero  and  assign 
a very  small  apriori  accuracy  estimate.  This  effectively  enforces 
the  desired  conditions.  This  approach  to  incorporating  splines 
in  a le^st  squares  adjustment  was  developed  by  DBA  under  a contract 
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for  ETL  on  the  ability  to  perform  triangulations  with  Side 
Looking  Radar  (DBA  1974)  . 

The  preceding  paragraohs  outline  the  method  for 
introducing  spline  constraints  into  a photoqrammetric  adiustment. 

The  following  paragraphs  will  use  OPS  data  as  an  example  of  how 
these  constraints  could  be  useful. 

As  was  discussed  earlier,  the  OPS  receiver  can  be 
programmed  to  output  positional  data  at  any  desired  time.  In 
Section  4. 6. 1.3  on  Error  Modeling,  it  was  assumed  that  the  receiver 
was  triggered  by  an  interrupt  which  was  generated  by  the  camera 
system  to  yield  position  at  the  time  of  the  exposure. 

An  alternative  to  this  ,is  to  record  the  OPS  position 
data  at  the  highest  possible  frequencv.  This  recording  would  be 
accomplished  completely  independent  of  the  photonraohv  process. The  only 
requirement  being  that  the  time  of  each  exposure  also  be  recorded. 

» 

The  GPS  positional  observations  would  be  input  directly 
into  the  photogrammetric  adiustment,  where  they  would  affect  only 
the  recovery  of  the  spline  parameters.  Since  these  spline  param- 
eters are  included  in  the  revised  colinearity  equations  (Eq . 4 . 6 . 2-5)  , 
the  overall  affect  is  to  constrain  the  photographs  to  lie  on  the 
flight  path  defined  by  the  continuous  GPS  data.  If  there  are  any 
systematic  discrepancies  between  the  GPS  data  and  the  photogrammetric 
data,  they  can  be  compensated  for  by  the  use  of  the  standard 
error  modeling  techniques  described  in  Section  4.6.1. 

It  is  interesting  to  note  that  if  one  assumes  equal 
length  spline  intervals  and  a uniform  GPS  sample  rate,  some 
computational  savings  can  be  made.  This  results  from  the  fact 
that  the  normal  equation  contributions  for  the  GPS  are 


1 
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strictly  a function  of  time,  relative  to  the  midpoint  of  the  spline. 
With  the  above  assumptions,  it  is  apparent  that  the  normal  equations 
for  this  set  of  parameters  will  be  constant  and  as  such  can  be 
computed  prior  to  the  adjustment.  This  precomputed  matrix  can  be 
added  to  the  photogrammetric  contribution  of  the  normal  equations 
during  the  formation  algorith.  This  is  not  true  of  the  correction 
vector  (c)  however.  The  GPS  contribution  must  be  computed  for  this 
vector  because  it  is  a function  of  time  and  the  observational  dis- 
crepancy vector.  This  principal  of  precomputing  normal  equations 
has  been  demonstrated  by  DBA  1976,  in  the  Offline  Orthoprinter 
Software  (OOPS)  development.  This  principal  has  also  been  employed 
by  geodesists  for  years  in  least  squares  reductions.  Geodesists 
compute  normal  equations  once  and  reuse  them  for  subsequent  iterations 
updating  only  the  constant  column  each  time.  In  OOPS,  a set  of 
normal  equations  were  set  uo  and  used  for  other  data  sets  and 
recomputed  only  if  the  result  indicated  it  was  required. 

With  the  technique  outlined  above, one  would  precompute 
only  a portion  of  the  normal  equations  and  add  it  to  the  rigorously 
formed  photogrammetric  normals. 

Earlier  in  this  section  an  assumption  was  made  that  w~ 
would  consider  only  a single  strip  of  photography.  The  immediate 
question  arises, "why  not  a block  of  photoqraphy?"  There  are  two 
basic  approaches  to  using  splines  for  a block  of  photoqraphs. 

The  first  alternative  is  the  reinitialization  of  the  spline  function 
for  each  strip,  i.e.  do  not  enforce  the  continuity  constraints 
between  the  last  spline  of  strip  1 and  the  first  spline  of  strip  2, 
and  do  not  use  the  GPS  data  for  the  time  elapsed  between  strip  1 
and  2 . 
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The  second  alternative  is  to  carry  the  spline  functions 
around  from  one  strip  to  the  next.  Obviously,  there  will  be  no 
photogrammetric  data  available  between  the  end  of  strip  1 and 
the  beginning  of  strip  2,  however,  the  continuous  GPS  observations 
should  be  sufficient  to  derive  the  required  parameters. 

This  discussion  of  the  use  of  splines  to  model  GPS 
data  is  only  one  possible  application  of  splines  in  a photogrammetric 
type  reduction.  Two  additional  areas  which  could  effectively  use 
this  type  of  modeling  are;  (1)  the  use  of  inertial  surveying  control 
information  and  (2)  in  conjunction  with  the  geometric  constraints 
described  in  Section  4. 6. 2. 4. 
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Banded-Bordered  Structure 
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The  banded-bordered  structure  of  the  photogrammetric 
normal  equation  coefficient  matrix  is  a developmental  outgrowth 
of  the  strictly  banded  coefficient  matrix.  This  development  of 
a banded-bordered  system  was  necessarv  to  accommodate  various 
classes  of  parameters  previously  not  considered  or  parameters 
which  would  detract  from  the  purely  banded  structure  of  the  N or 
N portions  of  the  coefficient  matrix.  Classes  of  parameters  are 
generally  considered  as  those  parameters  relating  to  a specific 
group  of  observations.  For  example:  GPS  Data,  as  previously 

descibed,  may  be  common  to  all  photos  in  the  block  (Error  Model 
Parameters,  Section  4. 6. 1.3)  or  may  be  common  to  portions  of 
strips  or  the  entire  strips  (Spline  Formation,  Section  4. 6. 2. 5). 

Miscellaneous  parameters  pertaining  to  certain  sets  of  observations 
(i.e.  Equal  Elevation  Constraints,  Section  4. 6. 2. 3)  may  also  be 
placed  in  the  border.  Such  information  would  normally  detract 
from  the  banded  structure  of  the  N matrix  by  forminq  off-diagonal 
terms.  With  the  subsequent  development  of  recursive  partitioning 
for  solving  banded-bordered  systems,  the  maintenance  of  a minimum 
bandwidth  is  essential.  Thus,  it  may  be  realized  that  a border 
will  permit  the  handling  of  more  varied  data  than  previously  used 
for  the  normal  equation  solution. 

Throughout  this  section,  examples  of  the  banded-bordered 
structure  will  be  presented  utilizing  algorithms  previously  developed 
for  GPS  Auxiliary  Data  (Sections  4.6.1  and  4.6.2).  Subsection  4.6. 3.1 
will  consider  the  GPS  Error  Modeling  information  as  auxiliary  data, 
whereas.  Sections  4. 6. 3. 2 (Structured  Border)  and  4. 6. 3. 3 (Double 

Fold  Solution  Algorithms)  will  consider  the  GPS  spline  information 
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as  auxiliary  data. 


4.6. 3.1  General  Formulation  of  Banded-Bordered  Normal  Equations 
In  the  formulation  of  a banded-bordered  system,  the 
initial  geometric  relationships  between  parameters  are  represented 
by  analytical  projection  of  arbitrary  three-dimensional,  object 
space  coordinate  system  into  a two-dimensional  image  space  system. 
Considering  the  GPS  Error  Modeling  Algorithms  discussed  in  Section 
4. 6. 1.3,  the  final  projective  equations  expressed  in  matrix  form 


+ B 6+6  8 + B 5 = e 


(4.6. 3-1) 


where; 


B 6 - Photo  parameter  coefficients  and  correction 

6 2 - Error  model  coefficient  and  corrections 


B 6 - Object  point  coefficients  and  corrections 


From  equation  (4. 6. 3-1)  the  following  total  set  of  observational  equations 


nay  be  derived; 


+ B6  + 6 8 + B6  = 


v - 6 


- 6 = e 


which  may  be  rewritten  as ; 


where  ; 


v + B 6 = e 


(4.6. 3-2) 


v ; B 


B ; 6 = 6 I ; e = 


-I  0 0 


0 -I  0 


0 0 -I 


('’.6.3-3) 
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In  addition  a diagonal  weight  matrix  is  formed; 


Brown  (1955)  has  shewn  that  the  least  squares  solution  for 
the  minimization  of  the  variance  is  obtained  when  the  vector  6 is  defined 
by  a system  of  normal  equations  of  the  form; 

6 = N-1  c (4. 6. 3-4) 

in  which; 

N = BT  w B (4.6. 3-5) 

-T 

c = B we 


Substituting  equation  (4.6. 3-3)  into  equations  (4.6. 3-5)  we  obtain; 
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By  redefinition  and  substitutions  into  equation  (4. 6. 3-4),  the  normal 
equations  become: 


<5 
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II 
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-T 
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ft 
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N 

N 

NJ 


(4. 6. 3-7) 


inhere  in: 

N - consists  of  the  relationships  between  N and  N, 

N - consists  of  the  relationships  between  ft  and  N,  and 

N - consists  of  the  relationships  between  N and  ft. 

The  N,  ft,  and  NT  terms  constitute  the  border  structure.  An  attempt  must  now 

be  made  to  solve  this  sizable  system  of  equations. In  an  effort  to  reduce  the 
inversion  effort  of  the  large  coefficient  matrix  a foldinq  operation  will  be 
performed.  This  svstem  of  normal  equations  will  first  be  partitioned  and 
redefined  as: 

6’|=  | N’  JN-  l"1  f c'  1 

(4. 6. 3-8) 
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We  shall  now  define  the  inverse  of  the  normal  equations  to  be; 

-1 


N_1  = M’  = 


M*  j N1 
_T'l  ~~ 
N j N 


ri'  M' 

_ml  ■■ 

M M' 


(4.6. 3-9) 


Since  M'is  the  inverse  of  N by  definition,  we  can  say: 
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(''.6.3-10) 


by  conforming  the  required  multiplication  of  the  left  hand  side  of  equation 
(4.6.3-10)  we  obtain: 


(4.6.3-11) 


N' 

M- 

+ 

N' 

-'T 

M 

= I 

N* 

M' 

+ 

N' 

M' 

= 0 

-T' 

N 

M' 

+ 

N 

-'T 

M 

= 0 

-T 

N 

M' 

+ 

N 

M' 

= I 

• -t 

from  which  we  can  solve  for  M,  M , M and  obtain: 

. . -'-1  ••  -'t  -1 

M’  = (N*  - N MS  ) 

M'  = - M'  N'  N_1 

M’  = N_1  + N_1  if  M'  N N"1 


(4.6.3-12) 

(4.6.3-13) 

(4.6.3-10 


We  may  now  define  auxiliary  form  in  an  effort  to  condense  this  formulation: 

(4.6.3-15) 


Let  us  define:  Q - N N 


r*  = n'  o' 

D/  nrbctituting  equations  (4.6.3-12)throucfi  (4.6.3-14)  into  equation  (4.6. 3-8) 
rnd  solving, we  obtain: 

S'  = M’  (c'  - N'  N_1  c)  (4.6. 3- 16a) 

. ••  --1  -'t  • 

S =M  c-N  M 6 
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or  with  the  use  of  auxiliaries,. 

S'  = M'  (c1  - Q'  c) 

6 = M-1  c - Q - S' 


(4.6.3-16b) 


It  should  be  noted  that  the  largest  inverse  now  required  is  the  size  of 
the  N'  portion  of  the  coefficient  matrix  only  and  not  the  full  coefficient 
matrix. 

Figure  4.6. 3-1  illustrates  the  structure  of  a typical  set 
of  normal  equations  generated  during  a photograirmetric  block  adjustment 

with  error  model  parameters.  Notice  that  in  this  exanple  the  N sulmatrix 

o _ 

is  highly  patterned,  however,  the  N,  N and  N are  completely  filled  in. 

The  next  Figure  (4.6. 3-2)  illustrates  the  structure  which  results  from  the 

. _i  . — 1 ••  ■■  -i  _t 

first  folding  operation.  This  matrix  is  the  M'  or  N'  -N  (N+W)  N 
which  must  be  inverted.  It  should  be  noted  here  that  the  border  is  con- 
sidered completely  filled  in.  This  structure  is  typical  of  the  reduced 
normal  equations  in  the  orbitally  constrained  reductions  LOSAT  and  PLODS. 
Haag  and  Hodge,  1972  describes  a modification  to  the  standard  recursive 

partitioning  algorithm  to  accomodate  the  border.  This  algorithm  requires 

o ~ -1  ~T 

that  the  lower  right  block  (N  - N (N+W)  N ) be  inverted  and  kept  in  core 
for  the  solution  process.  This  restriction  can  become  impractical  when 
one  considers  the  overall  potential  use  of  the  border.  The  next  section 
outlines  an  efficient  algorithm  for  solving  banded-bordered  structures 
which  have  sane  specialized  characteristics. 
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Normal  Equations  for  a typical  Four  Photo  Strip 

(only  upper  triangulation  portion  is  illustrated) 


4. 6. 3. 2 


Structured  Border 


In  order  to  decrease  the  size  of  the  N portion  of  the 
coefficient  matrix  for  conventional  inversion  and  solution,  or, 
to  decrease  the  number  of  elements  to  be  stored  in  the  border 
for  application  of  recursive  oartitioning , a method  of  double- 
folding has  been  developed.  The  algorithm  for  this  double- 
folding are  presented  in  Section  4. 6. 3. 3.  Before  double-folding 
can  be  efficiently  utilized,  the  border  of  the  coefficient  matrix 

must  be  structured  similar  to  the  N and  N portions.  In  other 

o 

words,  an  attempt  to  structure  the  N portion  in  a symmetric  block 
diagonal  fashion  must  be  made  so  that  off-diagonal  elements  are 
eliminated.  It  is  evident  that  only  certain  classes  of  param- 
eters lend  themselves  to  such  structuring.  In  an  effort  to 
attain  a block  diagonal  formation,  the  auxiliary  data  must  be 
analyzed  and  the  elements  of  the  N portion  if  necessarv  to 
achieve  this  formation.  Examples  of  such  auxiliary  data  are  the 
GPS  spline  algorithms  developed  in  Section  4. 6. 2. 5.  Utilizing 
these  algorithms  and  following  the  developmental  format  presented 
in  Equations  (4.6. 3-4)  through  (4.6. 3-7)  the  system  of  normal  equa- 
tions in  the  form  of  equation  (''.6.3-7)  is  arrived  at.  This 
system  can  be  double  partitioned  as  follows: 


An  example  of  this  double  partitioned  coefficient  matrix  with  a 

structured  border  is  shown  in  Figure  4. 6. 3-3.  In  addition  to 

o 

decreasing  the  storage  required  for  the  N portion,  one  also 
reduces  the  computational  effort  required  since  only  the  block 
diagonal  form  need  be  computed. 

4. 6. 3. 3 Double  Fold  Solution  Algorithms 

With  respect  to  Equation  (4. 6. 3. 8),  the  single  folded 
algorithm  effectively  places  the  influence  of  the  N and  N'  portions 

of  the  coefficient  matrix  into  the  N'  portion.  Double  folding 

o 

will  not  place  the  effects  of  the  once  folded  N portion  of 
Equation  4. 6. 3-7  into  the  N portion. 

Again  the  following  definitions  are  made  which  refer 
to  Equation  4. 6. 3-7; 
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From  equations  4. 6. 3-9  and  4.6. 3-8  is  derived; 


N 

Z » 

I 

-l 

•z 

M 

'T 

Ln 

ft. 

II 

ft. 

(4.6.3-17) 


(4.6.3-1G) 


From  equations  4. 6. 3-7  and  4.6.3-12  may  be  related  the  elements  of  equation 
4.6.3-18  such  that; 


M'  = 


• _ 

N - N N N 

~m  ~ _ I _rn 

N1-  N N N 


_ =T  -1 
N - N N N 

ft  - N N_1NT 


(4.6.3-19) 
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By  relating  the  elements  in  equation  (4.6.3-19)  to  equation 
(4. 6. 3-9)  and  by  substituting  into  equations  (4.6.3-12)  thru 
(4.6.3-14)  the  following  is  obtained; 


M = (N-N  N-1  NT)  - (N  - N N 1 NT)  (ft-N  N_1  NT)  (NT-ft  N 1NT)[  _1 


. ' - =T  • = •— 1 -1 

M = - M (N  - N N N ) (N  - N N N ) 

ft  = (S  - N N-1  NT)-1  + (ft  - N N-1  NT)_1(NT-S  N-1NT) 

. ~ "-1  =T  n = "-1  -T  -1 

M(N-N  N N ) (ft  - N N N ) 


(4.6.3-20) 


The  remaining  elements  of  the  inverted  coefficient  matrix  may  now 
be  defined  by; 


N'  = N 


(4.6.3-21) 


and  by  substituting  equation  4.6.3-21  and  4.6.3-18  into  equation 
4.6.3-13 


M'=-[m  Ml  [n]  N-1  = (yP  + M N)  N_1  =|M  (4.6.3-22) 

I ft  I IN  “ (MT  N + M N)  N_1  | ft 


therefore ; 


M = - (MNN1+MN  N_1) 
M = - (MT  ii  n’1  + I i N-1) 


(4.6.3-23) 


Similarly,  equations  4.6.3-21  and  4.6.3-18  can  be  substituted  into 
equation  4.6.3-14  to  obtain; 

M = N_1  + N-1  [nT  NT]  M M N N-1 

MT  ft  N_ 

"-1  "-1  -T  • - “-1  "-1  ~T  'T  - "-1  "-1  -T  ' = "-1 

= N + N N MNN  +N  N K NN  +N  N MNN  + 


N-1  NT  ft  ft  N_1 
-101- 


(4.6.3-24) 


I 
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Relating  elements  from  equations  4.6.3-20,  4.6.3-23  and  4.6.3-24 
to  corresponding  elements  in  equations  4.6.3-12  through  4.6.3-14 
a new  set  of  auxiliary  may  be  defined. 


"-1  -T 
Q = N N 


"-1  =T 
Q = N N 


r = nn1nt  = nq 


(4.6.3-25) 


= --1  =t  = ' 

R = N N N = N Q 


r = nn1nt  = nq 


o o - -1  ~T  =T. 
Q = (N  - R)  (N  - R ) 


R = (N  - ft)  (N  - R)-1  (N  - ft)T 


(N  - R)  Q 


By  substituting  into  equations  4.6.3-20,  4.6.3-23  and  4.6.3-24 
the  equations  can  be  simplified  to  be; 


. . o -1 

M = (N  - R - R) 

o ,o  ',-1  o . ot 
M = (N  - R)  + Q M Q 

. 0»T» 

M = - M Q 


- . T ~ 

M = - (M  Q +MQ  ) 

~ 'rp  t O 'T 

M = - (M  Q + M Q ) 


M = N’  1 + O (ft  - R)-1  QT  - (O  - O 6)  M (Q  - 0 5) 


(4.6.3-26) 
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The  result  of  this  second  foldina  is  that  the  largest 

. o 

matrix  now  required  for  conventional  inversion  is  the  N or  N 
portions  of  the  coefficient  matrix.  To  illustrate  the  potential 
utility  of  this  algorithm,  two  cases  of  GPS  data  with  the  spline 
modeling  described  in  Section  4. 6. 2. 5 shall  be  considered.  The 
first  case  will  have  the  splines  in  the  border  and  in  the  second 
case, the  photo  orientation  elements  will  be  in  the  border  with 
the  splines  in  the  N position. 

Figure  4.6. 3-3  illustrates  the  normal  equations 
generated  by  a typical  four  strip  block  of  photographs  where  each 
strip  contains  20  photographs.  The  time  span  for  each  spline  is 
arbitrarily  defined  as  five  photographs  in  length.  The  N or  N 
matrices  are  not  completely  illustrated  but  they  do  represent 
the  worst  case  condition. 
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The  N submatrices  represent  the  photo  orientation 

angles  arranged  by  the  standard  cross  strip  numbering  scheme. 

o 

N is  the  submatrix  generated  by  the  spline  coefficients.  There 
are  basically  three  elements  to  this  matrix  which  are  illustrated 
by  the  numbers  1,  2,  and  3 in  the  figure.  The  first  is  the  con- 
tribution to  the  spline  parameters  by  the  GPS  data  (1)  in  the 
figure.  Elements  2 (in  the  figure)  are  the  contribution  to  the 
spline  coefficients  from  the  photogrammetr ic  exterior  orientation 
parameters  and  elements.  The  third  element  (in  the  figure) 
represents  the  contribution  from  the  spline  continuity  constraints. 
Elements  1 and  3 are  those  which  were  previously  discussed  for 
potential  precomputing  of  their  values. 

It  is  important  to  notice  the  hiqhly  patterned  structure 
of  the  N,  N and  N matrices  which  is  in  contrast  to  the  discussion 
in  Section  4. 6. 3-1. 

Figure  4. 6. 3-4  illustrates  the  results  of  tne  first 

. o - 

folding  operation  when,  the  N portion  is  folded  into  the  N,  N N 

submatrices.  Notice  again  that  inthis  figure,  N is  still  highly 

o 

patterned  and  N has  a banded  structure. 

0 

The  second  folding  (N  into  N)  requires  the  inverse  of 

o - 

the  matrices  (N-R)  illustrated  as  the  lower  riqht  block  of  Figure 
4.6. 3-4.  Naturally,  this  inversion  process  will  cause  the  matrix 
to  be  completely  filled,  however,  if  one  can  assume  that  the 
elements  beyond  the  bandwidth  are  insignificant  numerically, 

then  the  inverse  will  have  the  same  essential  structure  as  the 

* 

normal  equations  . Based  on  this  assumotion, the  results  of  this 
* This  assumption  may  well  be  invalid  and  requires  further  testing. 
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FIGURE  4. 6. 3-3 

Typical  Normal  Equations  for  a 4x20  Photo  Block  with 
Five  Photo  Splines 


FIGURE  4. 6. 3-4 

Results  of  the  First  Fold  on  the  Normal  Equations 
(Case  I) 


second  folding  are  illustrated  in  Fiqure  4. 6. 1-5.  Note  that  in 
this  figure  the  matrix  is  also  a banded  structure  with  a very  large 
bandwidth  compared  to  the  original  (Figure  4.6.3-/1).  However 
recursive  partitioning  can  be  emploved  again  to  solve  this  system 
of  equations. 

Figure  4.6.3-6a  represents  the  normal  equations 

resulting  from  the  Case  II  examole,  namely,  the  interchange  of 

. o 

the  N and  N portions  of  the  normal  equations.  The  same  double 
folding  algorithm  is  again  employed  to  solve  this  system  of 
equations.  The  results  of  this  operation  are  illustrated  in  Fiqure 
4. 6. 3-3b. 

There  are  manv  additional  wavs  to  arrange  the  data 

o 

within  the  N and  N matrices  which  effect  the  structure  of  the 
first  folding  and  second  foldina.  Much  more  analvsis  is  needed 
to  determine  the  utility  of  this  procedure  in  practice. 
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4.6.4 


Autoregressive  Modeling 


In  Section  4.6.1  a detailed  description  of  how 
various  systematic  error  models  could  be  directly  introduced 
into  the  photogrammetric  block  adjustment  formulation  were 
presented.  All  of  these  models  assumed  that  the  basic  observations 
were  independent  from  one  another  or  were  functionally  related. 
There  is  another  class  of  observations  which  are  occasionally 
available  for  the  photogrammetric  problem.  These  are  observations 
which  are  serially  correlated.  The  technique  employed  in  SAGA, 
(Brown  & Trotter,  1969)  called  Autoregressive  Feedback  for 
including  serially  correlated  range  observations  can  also  be 
employed  in  the  photogrammetric  block  adjustment. 

This  technique  is  particularly  well  suited  for 
observations  which  are  time  dependent  and  occur  at  a very  high 
frequency  such  as  the  GPS  positional  data  discussed  in  Section  4.3.1 
or  the  inertial  survey  data  described  in  Section  4.4.2. 

The  following  sections  present  the  basic  mathematical 
development  of  the  autoregressive  feedback  technique  with  emphasis 
on  the  GPS  type  observations. 

4 . 6 . 4 . i The  Autoregressive  Model 

A stationary  sequence  of  serially  correlated  errors 
c i , r^-^,  2'  • • ■ is  qoverned  by  an  autoreqressive  process  in 

which 

el  = V- 1-1  + "2'  1-2  + •••  + ap  rl-p  + nl 

4 . 6 . 4 - 1 

where  «'s  are  constant  coefficients  and  is  a random  impulse 

2 

of  zero  mean  and  variance  o . According  to  this  process,  the 

i th  error  in  the  sequence  is  generated  as  a fixed  linear 
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1 


combination  with  a superimposed,  strictly  random  impulse.  The 
process  indicated  in  equation  (4. 6. 4-1)  is  said  to  be  the  of  pth 
order  because  p coefficients  are  involved  in  its  description. 
Specific  examples  of  the  character  of  errors  generated  by  various 
first  order  autoregressive  processes  are  to  be  found  in  Brown,  Bush, 
Sibol  (1963).  DBA's  past  experience  with  autoregressive  feedback 
on  systems  similar  to  GPS  indicates  that  a low  order  process  is 
probably  sufficient  to  model  the  serial  correlation  in  GPS  data. 

4. 6. 4. 2 Inverse  Covariance  Matrix  of  Autoregressive  Process 
As  pointed  out  in  earlier  work,  the  practical 
utility  of  the  autoregressive  process  in  the  adjustment  of 
observations  stems  from  the  fact  that,  for  a given  set  of  auto- 
regressive coefficients  (which,  can  be  estimated  from  the 
observations) , the  inverse  of  the  covariance  matrix  of  the 
observations  can  be  computed  analytically.  Equally  important, 
the  inverse  is  a multidiaqonal  matrix,  in  which 
the  number  of  diagonals  are  equal  to  2p+l  for  a process  of  the 
p-th  order.  In  Brown,  Bush,  Sibol  (1964),  it  was  demonstrated 
that  the  basic  result  of  this  effect  (derived  oricrinally  by 
Wise,  1955)  could  be  put  into  a more  convenient  form.  Namely, 
if  A denotes  the  covariance  matrix  of  an  autoregressive  process 
governed  by  coefficients  a . , .*  7 , . . . , , then  A ^ can  be  expressed 

as  the  followinq  product  of  lower  and  upper  trianqulatar  matrices; 
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in  which  an=0  for  n>p.  If  this  result  is  applied,  for  instance 
to  a third  order  process  involving  coefficients  a]/a2,a3'  one 
finds,  upon  performing  the  above  matrix  multiplication  that  A ^ 
assumes  the  form 


o A = 


a_.'(  b_~  c.j  d \0  0 0 0 ...  0 0 0 0 

bus  a.,  b.;  c dN  0 0 0 ...  0 0 0 0 

\ 


c_i  b_i  a_i  bed  >^0  0 


0 0 0 0 


d c b a b c d\0...0  0 0 0 

O'^d  c b a b c d\. . . 0 0 0 0 


0 0 \d  c b a b c*.\0  0 0 0 

0 0 0v^cbab...0'v00  0 

0 0 0 0 Nd  c b a ...  0 0 \0  0 

• • .V  . 

\ • • • • • < . 

••••••••  I • • 

00000  6^  00. ..abed 

0 0 0 0 0 0 0 ^0  ...  b o-i  b-i  c-i 

0 0 0 0 0 0 0 0\...  C b_i  a_o  b_rl 

00000000  C-i  I)..,  a.. 


ooooo  6^  oo 
0000000*  0 


in  which; 


i 

! 


a ~ 1 

+ cf,7  + <y27  + 

°-i 

= 1 + cr,7  + Ct] 

°-2 

= 1 + o',2 

a-3  = 1 

b = -O', 

+ <7,  a2  + a2a3 

b-, 

= “O',  + O',  fV2 

b-2 

= ~al 

C - -Ct'2 

d 

+ <v2o3 

C-l 

- -a 
2 

4.6 .4-4. 


This  illustration  is  sufficiently  general  to  demonstrate  that 
certain  properties  of  the  inverse  are  of  pivotal  importance  to 
the  concept  of  autoregressive  feedbacks.  The  properties  are; 

(1)  the  inverse  is  multi-diagonal,  the  number  of 
diagonals  being  equal  to  7 for  p=3  (or  equal 
to  2p+l  in  general). 

(2)  except  for  pxp  submatrices  comprising  the  upper 
and  lower  corners  of  the  matrix,  the  elements 
of  each  diagonal  are  constant,  the  number  of 
different  constants  being  equal  to  p,  the  order 
of  the  process. 

(3)  the  formulas  for  the  elements  of  the  matrix  are 
subject  to  a systematic,  easily  generalized 
development . 

Expressions  for  second  and  first  order  processes  can  be  obtained 
from  the  above  development  by  sucessively  equating  oc^  an<3  a2 
equal  to  zero  in  equation  4. 6. 4-4 

4. 6. 4. 3 Estimation  of  Autoregressive  Coefficients 

Let  us  assume  that, by  some  means  one  has  obtained 

a vector  of  residuals  (v. , v_,...,v  ) that  constitutes  an 

1 z n 

estimate  of  the  vector  of  actual  errors  (f,,c_,...,r  ).  One 

12  n 

can  then  generate  a set  of  autocorrelation  coefficients  p^,p0... 
in  the  usual  manner  from; 


E vi  vi-* 

P*  =■ 


4 .6 . 4-5 
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V!  =0|  Vl-1  + <*3  vt_3  % + rj  j 


and  this  can  be  regarded'  as  defining  a set  of  observational  equations 


! 

involving  the  autoregressive  coefficients  as  parameters.  The 
least  squares  adjustment  technique  can  be  employed  to  generate 
a system  of  normal  equations  for  the  determination  of  the  oi's. 

These  turn  out  to  assume  the  form; 


[-  — 1 
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— 
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In  practice,  it  is  desirable  to  establish  the  autoregressive 
process  of  lowest  order  which  satisfactorily  models  the  observed 
process.  This  can  be  accomplished  by  beginning  with  a first 
order  process  which,  from  equation  4. 6. 4-7  would  lead  to  the 
estimate  ; 


<v,  = Pi  . 


4. 6. 4-8 
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; 


; 

: 


i 


1 


Using  this  estimate  of  the  secondary  residual  n.  can  be 
computed  as: 

= *1  - O',  Ci-  | . 4. 6. 4-9 

If  these  secondary  residuals  turn  out  to  be  serially  uncorrelated 
the  first  order  process  can  be  accepted  as  satisfactory.  To 
determine  whether  or  not  the  n's  are  independent,  their 

or<3er  coefficient  of  correlation  will  be  computed  from; 

_ V l-j 

r~ . 4.6.4-10 

Zv* 

The  statistical  signifance  of  r can  be  tested  using  the  well 
known  result  that  the  quantity  ; 


z = £ 


Ur 

1-r 


4.6.4-11 


is  approximately  normally  distributed  with  mean  and  variance  of 


Z 


= i 

c 


In 


Up 

UP 


4.6. 4-12 


= V(n-3)  4.6.4-13 

in  which  p is  the  true,  but  unknown  correlation  coefficient. 

It  follows  that,  at  the  95%  level  of  confidence,  r will  differ 
significantly  from  zero  only  if 


2 > 


^ • ^6/  Jn-3  t 


4.6.4-14 
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If  r should  turn  out  to  be  significantly  different  from  zero. 


then  an  attempt  to  model  the  system  using  a second  order  process 
would  be  initiated.  Here,  according  to  equation  4. 6. 4-7,  the 
coefficients  would  be  established  from 


1 P, 

P,  1 


4.6.4-15 


which  has  the  solution 


a,  = (Pi  - Pi  pa)A]-r%)> 

os  = (-p?  +Pa)/0-*$)- 


4.6.4-16 


The  secondary  residuals  for  the  second  order  process  are 


<r,  - 


4.6.4-17 


The  serial  independence  of  these  residuals  would  be  tested  in 
precisely  the  same  manner  as  described  above  for  the  first  order 


process . 


The  above  procedure  would  be  repeated  until  a 


satisfactory  fit  is  obtained.  Once  a satisfactory  autoregressive 
model  has  been  established,  the  spectral  density  function  can 
be  computed  analytically  from  the  following  elegant  result 


derived  from  Wise: 


1/(0)  = 0?/(- 1+  or,e,0+  <v2e‘20+  ...  + a ’ e‘P°)  H + 0',o",G+  "2e  ’2°+...  + 0'  c 'r0) 


4.6.4-18 


where  i = /-I  and  0 assumes  the  discrete  sample  values, 


2rr  6ir 


TT,  27T  . 


4.6.4-19 
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If  At  denotes  the  time  between  successive  points  and  T=nfit  denotes 
the  total  time  span,  9 may  be  put  in  the  form, 


0 


1 


2 lAt 
~T~  ' 


i 1/  2, , * , ,n 


4.6.4-20 


where  6^  may  be  said  to  correspond  to  the  frequency  of  1/2^ at 
cycles  per  second.  For  the  first  order  process,  equation  4.6.4-18 
reduces  to  the  following  well  known  expression  for  the  spectral 
density  function  of  the  damped  exponential  autocorrelation 
function. 


y(0)  = . o */  (1  - 2 ° cos  G + D?) 


4.6.4-21 


4. 6. 4. 4 Refined  Normal  Equations 

From  the  foregoing, it  can  be  seen  that  a successful 
autoregressive  analysis  of  residuals  provides  the  solution  to 
two  central  problems  of  random  error  analysis:  (1)  the  determination 

of  the  inverse  covariance  matrix  of  the  observational  vector  and 
(2)  the  determination  of  the  spectral  density  function  of  the 
error  process.  It  remains  to  be  shown  precisely  how  this  information 
can  be  used  in  a refinement  of  the  adjustment  and  what  its  use 
entails  in  the  way  of  additional  complication.  For  this  purpose 
the  specific  case  of  an  adjustment  of  a single  channel  of  obser- 
vations qoverned  by  a first  order  autoreqressive  process  shall  be 
considered  . This  is  sufficient  to  demonstrate  the  general  principles 
of  the  operation.  Accordingly, 


Pf.-t  + *» 


4.6.4-22 
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define  a first  order  process  in  which; 


e (nt)  = 0 


4.6.4-23 


var  (V,  ) ~ (variance  of  high  frequency  component  of  noise).  4.6.4-24 


cov  (V^i-fc)  = 0 for  all  k>  0 


4.6.4-25 


Then  it  is  easily  shown  that; 

EU,)  = 0 > 

var  (^,)  ~ 0 ~ °2)  = total  error  variance. 


P*  O 
1 -O?  ' 


= P*. 


4.6.4-26 


4.6.4-27 


4.6.4-28 


4 .6.4-29 


It  follows  that  the  covariance  matrix  A of  an  n vector  of  errors 
is; 


A - 


1 -pr 


1 

P 

P2 
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P' 

P 
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, n-1 


>n-2 


4.6. 4-30 


pn  pn-1  pn-2 
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The  consequences  of  the  patterned  and  regular  structure  of  A 
shall  now  be  traced  through  the  least  squares  reduction.  First, 
note  that  A-1  can  be  decomposed  as  follows  . 


L. 


02A  - = (1  + PJ)I  - PU  - PU  - P2A 


4.6.4-34 


In  which  I is  an  nxn  unit  matrix  and  U and  A are  nxn  matrices 


of  the  form 


0 1 0 0 ...  0 
0 0 10  0 
0 0 0 1 0 


1 0 0 ...  0 

0 0 0 ...  0 

0 0 0 ...  0 
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0 0 0 0 ...  1 

0 0 0 0 ...  0 


0 0 0 ...  0 

0 0 0 ...  1 


N is  now  defined  as  the  coefficient  matrix  of  the  normal  equations, 

N = BV  B 

h 

= .1^  (bj  b]  ...  J)  [(1  + P?)  I - PU  - P UT  - P2  A ] b2 


4.6.4-36 


This  expands  to 


- * l(ltP’)(b'[  bk  + b,  +...+  b„  b„) 

CT  _ j 

-P  (b[  bK  bj  b,  +...+  b»-! 

-P  (bl  bk  H bj  b2  -*...+  bl  bB_j) 

-p3  (b|  bj  + b„  b0)| 


4.6.4-37 


If  we  define  , 


AL>!  = bt  (1-  b, 


4.6.4-38 


and  note  that  then  , 

bj_,  bj  = bj_,  (bj_j  + Abj.j)  = b,_j  b + b t x A^..,  , 


4.6. 4-39 


b|  bj..!  ~ bt  (b t - Ab j_ , ) - bj  b,  - bj  Ab,  , , 


4.6. 4-40 


it  follows  that  f 

bj-i  b j + bj  b j _i  = bj_j  bt_j  + bt  b,  4 (b t _ x - b,)  Ab,_j 


= b j _ j b(_j  4 bj  bj  + Abj_2  Abj_j 


4.6.4-41 


When  this  result  is  substituted  into  4.6.4-37  and  appropriate 
algebraic  manipulations  are  performed,  one  obtains 


N = — [ (1-2P4P*)CoJ'  b,  + bj  b8  4...+  bjbB) 

I 


4.6.4-42 


+ P (A bj  Ab,  + Abj  Ab2  -i . . Abj_2  Abn.,) 


+ p ( 1 - p) (b j b,  4 bj  bj 


which  is  the  same  as 


N 


= — { ( 1 - P) 1 BTB  4 0 ABT  ab  4 p(l-p)(bj  b,  4 bj  bn)j 


4.6.4-43 
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- - 


in  which 


ab=  : 


t>s  - bj 

^ ~ b2 


Abn-1  bc  - bn_j 


4.6.4-44 


In  an  analogous  manner  one  obtains  for  the  right  side 
of  the  normal  equations 

T -l 

c - B1  A e 

’ t 

= [ (l-ps)  BT  e + oABT  Ac  + p(l-p)(b[  Cj  + bB  eB)  J 

_3  l 


4.6.4-45 


in  which; 


Acn_! 


f3  “ 


Cn  ” CB_j 


4.6.4-46 


Equations  4.6.4-43  and  4.6.4-45  are  the  results 
sought.  In  analyzing  them, it  is  first  noted  that  when  >=0,  the 
normal  equations  reduce  to  the  usual  expression: 


— (B1  B)  6 = --  BT  c 
cr3  cra 


4.6. 4-47 


Next  , note  that  since  AB  and  Ar  are  of  first  order,  the 
expressions  AB  AB  and  AB  Ac  are  of  second  order  and  one  can 


assert  that  ; 
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btb  » abt  ab 


4.6.4-48 


» AB^  A € 

This  implies  that  as  long  as  p is  not  too  close  to  unity,  the 

normal  equations  will  be  dominated  by  the  leading  terms  and 

the  solution  vector  6 will  be  only  slightly  dependent  on  the 

value  of  p.  Hence,  in  some  situations  the  solution  is  but 

weakly  affected  by  the  presence  of  even  rather  moderate  serial 

correlation.  Be  this  as  it  may,  the  covariance  matrix  of  the 

solution  vector  is  very  strongly  dependent  on  the  degree  of 

T 

serial  correlation.  Then,when  B B dominates  the  normal  equations, 
the  covariance  matrix  of  the  parametric  vector  is  given  by 

(BTBf‘  . 4.6. 4-4 

O-0)3 

Thus  if  p were  actually  equal  to  0.9  and  one  were  to  ignore  this 

fact  in  the  adjustment,  the  solution  vector  itself  would  probably 

2 

not  be  very  much  affected,  for  the  factor  (1-p)  of  the  leading 

terms  on  both  sides  of  the  normal  equations  would  cancel  out  , 

but  the  covariance  matrix  associated  with  the  solution  would  be 

2 

incorrect  by  a factor  100  (or  l/(l-p)  ). 

Because  of  the  relative  insensitivity  of  the  solution 

vector  to  serial  correlation,  the  residuals  from  the  adjustment 

are  also  relatively  insensitive  to  serial  correlation  (remember 

an  adjustment  restricted  to  a single  observational  vector  r is 

considered  herein) . This  means  that  it  is  a sound  and  defensible 

practice  to  estimate  the  autoreqressive  function  from  residuals 
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of  a preliminary  adjustment  in  which  the  observational  errors 
are  initially  considered  to  be  uncorrelated.  The  autoregressive 
function  thus  initially  determined  can  then  be  fedback  into  the 
solution  according  to  the  development  of  this  section  and  the 
resulting,  more  nearly  correct  normal  equations  can  be  solved  to 
generate  fresh  residuals  for  a more  refined  autoregressive  analysis. 
Clearly,  this  process  can  be  iterated  to  stability  (ordinarily, 
a single  iteration  is  sufficient) . Thus,  autoregressive  feed- 
back is  an  adaptive  process  that  ultimately  becomes  independent 
of  initial  or  apriori  estimates  of  noise  variance. 

Returning  to  the  normal  equations  4.6.4-43,  4.6.4-45, 
it  is  noted  that  the  extra  computations  entailed  by  rigorous  con- 
sideration of  a first  order  autoregressive  process  are  surpris- 
ingly minor.  The  BTB  and  BTe  terms, being  of  similar  form,  follow 
the  same  logic.  Moreover,  as  is  clear  from  4.6.4-42,  the  equa- 
tions can  be  formed  in  a cumulative  manner:  the  i-1 st  step  of 

which  would  involve  only  observational  equations  from  the  i-th 
and  1 + 1 st  observations.  Hence,  the  formation  of  the  normal 
equations  for  a first  order  process  entails  a scheme  in  which 

a moving  pair  of  successive  observational  equations  are  processed 
at  each  step;  similarily,  a second  order  process  entails  a scheme 
in  \4iich  a moving  triplet  of  successive  observational  equations 
are  processed  at  each  step,  and  so  on.  A major  benefit  to  be 
derived  from  the  admission  of  an  autorenressive  process  into  the 
adjustment  is  the  attainment  of  more  realistic  results  from 

error  propagation.  Tn  particular,  it  would  permit  one  to 
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extract  whatever  advantages  are  to  be  gained  from  moderately 
high  sampling  rates  without  paying  the  usual  oenalty  of  an 
absurdly  optimistic  estimation  of  output  accuracies. 

The  process  outlined  here  is  very  general  in  nature 
and  can  be  applied  to  many  types  of  high  frequency  observations. 

DBA  has  successfully  implemented  this  algorithm  in  our  short 
arc  tracking  reduction  program, SAGA, for  handling  serially  correlated 
ranging  observations.  The  techniques  could  be  easily  adapted  to 
the  GPS  data  which,  if  collected  at  a high  frecruency,  could  exhibit 
high  serial  correlation. 


-125- 


4.6.5 


Automatic  Editin' 


As  pointed  out  in  Section  2.7  and  3.7  Automatic 
Editing  features  have  been  included  in  mnay  photogrammetr ic 
triangulation  programs.  However,  all  require  that  a new 
formation  and  solution  be  performed  each  time  one  or  more 
points  is  either  edited  or  reinserted. 

A technique  was  developed  by  Mikhail  and  Helmering  in 
1973,  whereby  points  could  be  edited  and  only  their  contribution 
subtracted  or  added  into  a set  of  previously  formed  normal 
equations.  Helmerina's  technique  has  been  exoer imentallv  tested 
on  the  TA-3/PA  at  DMAAC  and  was  demonstrated  to  save  a tremendous 
amount  of  time  in  the  real  time  editing  of  data.  However,  the 
technique  has  never,  to  our  knowledge,  been  applied  to  a ohoto- 
grammetric  block  adjustment.  Helmerinq’s  Doctorial  dissertation 
adequately  develops  the  mathematics  for  the  three  photo  alaorithm. 
This  development  is  easily  expanded  to  the  more  general  case 
and  will  not  be  presented  herein. 

Although  not  an  automatic  editing  scheme,  another 
technique  which  is  not  currently  utilized  in  photogrammetr ic 
triangulation  is,  interactive  editing  via  a CRT  display  screen. 

The  basic  concept  is  to  first  store  the  results  from  a trianqu- 
lation  on  either  a maqnetic  tape  or  disc  file.  This  tape  would 
then  be  input  to  an  editing  program .whereby  the  observational 
residuals  could  be  plotted  in  a varietv  of  formats  which  would 
be  operator  selectable.  The  ODorator  could  then  edit  a con- 
tribution of  this  observation  from  the  stored  normal  equation 

and  constraint  column  and  resolve  the  system.  The  results  would 
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then  be  displayed  for  further  editing. 

Another  automatic  editinq  technique  is  possible 
for  image  measurement  which  is  based  on  the  pair-wise  inter- 
section of  multiple  ray  points.  DBA  has  used  this  technique 
for  manual  data  editinq  very  successfully.  The  technique  shows 
tremendous  potential  in  automatically  detecting  image  measure- 
ment errors. 

The  combination  of  the  three  anproaches  outlined  above 
would  result  in  a orocedure  which  would  definitely  increase  the 
overall  throughput  of  a triangulation  process. 


5.0 


CONCLUSIONS  AND  RECOMMENDATIONS 


In  the  introduction  to  this  report  it  was  stated  that 
the  overall  objective  of  this  report  is  an  identification  of  new 
technology  which  could  either  increase  the  targeting  accuracy 
capabilities  for  the  MC&G  Community  or  to  increase  the  overall 
Photogrammetric  Triangulation  throughput.  A number  of  specific 
shortfalls  in  the  existing  triangulation  process  have  been  iden- 
tified, many  of  which  are  beyond  the  scope  of  this  study  to 
analyze  in  detail. 

The  primary  shortfall  identified  herein  is  the 
inability  of  current  techniques  to  exploit  all  of  the  information 
available  for  triangulation.  This  includes  auxiliary  sensor 
data  and  photographic  imagery  data.  In  general,  this  information 
can  be  included  into  an  adjustment  in  a number  of  methods;  direct 
observations,  functional  representation,  error  modeling  or 
statistically.  The  general  procedure  for  including  each  of 
these  types  of  models  are  included  herein  with  specific  examples 
of  what  types  of  information  could  be  included  with  each  model. 
Many  of  these  procedures  result  in  a set  of  normal  equations  which 
typically  exhibit  a banded-bordered  structure.  Although  solutions 
for  this  type  of  normal  equation  system  have  been  in  existence  for 
many  years,  nc  attempt  has  been  made  to  develop  an  efficient 
algorithm  which  exploits  the  specific  structure  of  the  border 
itself.  An  algorithm  for  this  is  developed  in  Section  4.6.3  of 
this  report. 

Since  the  objective  of  this  report  is  to  identify 

technologies  for  future  use, no  attempt  has  been  made  to  quantify 
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the  exact  benefits  which  one  could  expect  to  derive  by  implementing 


the  techniques.  Instead,  a set  of  alternatives  have  been  identified, 
many  of  which  need  much  more  analysis  and  testing  before  their 
specific  benefit  can  be  realistically  assessed.  Listed  below 
are  areas  which  DBA  believes  worthy  of  future  indepth  study. 

1.  Simulation  testing  of  the  various  techniques 
for  incorporating  GPS  positioning  data. 

a.  Error  Modeling 

b.  Spline  type  functional  modeling 

c.  Autoregressive  feedback  modeling 

2.  Implementation  and  testing  of  the  banded-bordered 
solution  algorithm. 

3.  Automatic  reordering  of  the  normal  equations 
with  emphasis  on  handling  the  border. 

4.  Geometric  constraints  bv  the  method  of  continuous 
traces . 

5.  Use  of  Mini-computers  and/or  interactive  systems 
for: 

a.  Material  availability 

b.  Coverage  Planning  and 

c.  Automatic  editing. 

Other  techniques  which  are  identified  in  this  report  can  be 
implemented  without  further  study.  These  include: 

1.  Distance,  Azimuth  and  elevation  constraints,  and 

2.  Camera  Calibration  and  film  deformation  models. 

A single  study  could  possibly  cover  the  first  three 

items.  In  that  during  the  GPS  spline  modeling  techninue  phase 

a banded-bordered  set  of  normal  equations  can  be  formulated  as 

in  Section  4.6.3.  This  system  could  then  be  solved  bv  use  of 

standard  methods  as  well  as  the  double  fold  technique.  The 

results  could  then  be  compared  with  respect  to  accuracy  and  run 
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times.  If  the  double  fold  techniaue  is  practical,  automatic 

reordering  algorithms  can  be  investigated  to  minimize  the  run 

c 

time  by  rearranging  the  N and  N elements  to  maximize  the 
efficiency  of  the  algorithm. 

Another  simulation  type  study  should  be  undertaken 
to  develop  and  test  the  algorithm  for  aeometric  constraint  type 
modeling  using  continuous  traces.  This  study  could  include  a 
small  set  of  three  photographs  of  real  data  to  validate  the 
approach. 

A third  area  relates  directly  to  the  overall 
triangulation  process  throughput.  This  is  the  use  of  interactive 
mini-computer  stations  for  analvzincr  the  Material  Available, 
Coverage  Planning  and  Editing  aspects.  As  pointed  out  in  Sections 
2 and  3,  these  steps  are  currently  performed  manually  and  are 
very  tedious  and  time  consuming.  With  the  current  technology 
and  the  low  cost  of  small  computers,  there  exists  the  potential 
to  greatly  reduce  the  tedium  and  time  through  automation. 
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METRIC  SYSTEM 


BASE  UNITS: 


Quantity 

Unit 

SI  Symbol 

Formula 

length 

metre 

m 

mass 

kilogram 

kg 

time 

second 

s 

electric  current 

ampere 

A 

thermodynamic  temperature 

kelvin 

K 

amount  of  substance 

mole 

mol 

luminous  intensity 

candela 

cd 

SUPPLEMENTARY  UNITS: 

plane  angle 

radian 

rad 

solid  angle 

steradian 

sr 

DERIVED  UNITS: 

Acceleration 

metre  per  second  squared 

m/s 

activity  (of  a radioactive  source) 

disintegration  per  second 

(disintegration)/a 

angular  acceleration 

radian  per  second  squared 

rad/s 

angular  velocity 

radian  per  second 

rad/s 

area 

square  metre 

m 

density 

kilogram  per  cubic  metre 

kg'm 

electric  capacitance 

farad 

F 

A-s/V 

electrical  conductance 

siemens 

S 

A/V 

electric  field  strength 

volt  per  metre 

V/m 

electric  inductance 

henry 

H 

V-s/A 

electric  potential  difference 

volt 

V 

W/A 

electric  resistance 

ohm 

V/A 

electromotive  force 

volt 

V 

W/A 

energy 

joule 

1 

N-m 

entropy 

joule  per  kelvin 

I'K 

force 

newton 

N 

kg-m/s 

frequency 

hertz 

Hz 

(cycle)/* 

illuminance 

lux 

lx 

im/m 

luminance 

candela  per  square  metre 

cd/m 

luminous  flux 

lumen 

Im 

cd-sr 

magnetic  field  strength 

ampere  per  metre 

A/m 

magnetic  flux 

webei 

Wb 

Vs 

magnetic  flux  density 

tesla 

T 

Wb/m 

magnetomotive  force 

ampere 

A 

power 

watt 

W 

J/» 

pressure 

pascal 

Pa 

Nm 

quantity  of  electricity 

coulomb 

C 

A-s 

quantity  of  heat 

joule 

1 

N-m 

radiant  intensity 

watt  per  steradian 

W/sr 

specific  heat 

joule  per  kilogram-kelvin 

)/kg-K 

stress 

pascal 

Pa 

N/m 

thermal  conductivity 

watt  per  metre-kelvin 

W/m-K 

velocity 

metre  per  second 

m/s 

viscosity,  dynamic 

pascal-second 

Pa-s 

viscosity,  kinematic 

square  metre  per  second 

m/s 

voltage 

volt 

V 

W/A 

volume 

cubic  metre 

m 

wavenumber 

reciprocal  metre 

(wave)/m 

work 

joule 

1 

N-m 

SI  PREFIXES: 

Multiplication  Factors 

Prefix 

SI  Symbol 

1 000  000  000  000  - 10'* 

(era 

T 

t 000  000  000  = 10* 

K'K* 

c; 

1 000  000  It)* 

mega 

M 

1 000  = 10’ 

kilo 

k 

100  =*  10a 

hecto* 

h 

10  = 10' 

deka* 

da 

0 1 ^ 10*' 

deci* 

d 

001  * 10"* 

c.entl* 

c 

0 001  10“  * 

mill! 

in 

0 000  001  * 10*  • 

micro 

0 000  000  001  10  * 

nano 

n 

0 000  000  000  001  10  '* 

plro 

P 

0 .000  000  (KK1{HK)  001  10  M 

femto 

1 

0 000  000  000  000  OOf)  001  to  •" 

afto 

a 

* To  bp  avoided  where  possible 


